{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# CS 237 Spring 2021, HW 07 \n", "\n", "#### Due date: Wednesday March 17th at Midnight (1 minute after 11:59pm on 2/11) via Gradescope (30 hour grace period)\n", "\n", " Late policy: You may submit the homework up to 48 hours late for a 10% penalty. Hence, the late deadline is Friday 3/19 at Midnight (with the normal 6 hour grace period). \n", "\n", "#### General Instructions\n", "\n", "Please complete this notebook by filling in solutions where indicated. \n", "\n", "For full credit, please take careful note of the following requirements:\n", "\n", "- Do NOT use any HTML tags in your notebook, as Gradescope will ignore them;\n", "\n", "- Do NOT answer questions by including images, as Gradescope will ignore them; and \n", "\n", "- You MUST \"Restart and Run All\" from the Kernel menu before submitting to Gradescope.\n", "\n", "**Any assignments which do not follow these requirements will not receive full credit.** \n", "\n", "\n", "\n", "There are 8 analytical problems and 4 programming problems. An introductory video will be posted on YT for\n", "the analytical problems, and the programming problems will be covered Friday in lab. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Here are some imports which will be used in code that we write for CS 237\n", " \n", "\n", "# Imports potentially used for this lab\n", "\n", "\n", "import matplotlib.pyplot as plt # normal plotting\n", "import numpy as np\n", "\n", "from math import log, pi,log,floor # import whatever you want from math\n", "from random import seed, random\n", "from scipy.special import comb\n", "from collections import Counter\n", "\n", "%matplotlib inline\n", "\n", "# Calculating permutations and combinations efficiently\n", "\n", "def P(N,K):\n", " res = 1\n", " for i in range(K):\n", " res *= N\n", " N = N - 1\n", " return res\n", " \n", "def C(N,K): \n", " return comb(N,K,True) # just a wrapper around the scipy function\n", "\n", "\n", "# Useful code \n", "\n", "def show_distribution(outcomes, title='Probability Distribution'):\n", " num_trials = len(outcomes)\n", " X = range( int(min(outcomes)), int(max(outcomes))+1 )\n", " freqs = Counter(outcomes)\n", " Y = [freqs[i]/num_trials for i in X]\n", " plt.bar(X,Y,width=1.0,edgecolor='black')\n", " if (X[-1] - X[0] < 30):\n", " ticks = range(X[0],X[-1]+1)\n", " plt.xticks(ticks, ticks) \n", " plt.xlabel(\"Outcomes\")\n", " plt.ylabel(\"Probability\")\n", " plt.title(title)\n", " plt.show()\n", "\n", "# This function takes the range and PMF of a discrete RV and draws the distribution. \n", "\n", "def draw_distribution(Rx, fx, title='Probability Distribution for X'):\n", " plt.bar(Rx,fx,width=1.0,edgecolor='black')\n", " plt.ylabel(\"Probability\")\n", " plt.xlabel(\"Outcomes\")\n", " if (Rx[-1] - Rx[0] < 30):\n", " ticks = range(Rx[0],Rx[-1]+1)\n", " plt.xticks(ticks, ticks) \n", " plt.title(title)\n", " plt.show()\n", " \n", "def round4(x):\n", " return round(x+0.00000000001,4)\n", "\n", "def round4_list(L):\n", " return [ round4(x) for x in L]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analytical Problem Instructions\n", "\n", "The Problems 2, 3, and 4 ask you to \"describe\" a random variable, which means:\n", "\n", "> (i) Give $R_X$ (you may schematize it if it is very complicated or infinite);
\n", "> (ii) List out the values of $f_X$ corresponding to each element of $R_X$;
\n", "> (iii) Draw a probability distribution, using the function draw_distribution provided in the previous cell.
\n", "\n", "As always, round to 4 decimal places at the last stage, using the functions round4(...) and round4_list(...) given above.\n", "\n", "A nice way to approach these is to do any complicated calculations in Python and then if you have\n", "to change something you won't have to redo all the calculations. Plus, you will make fewer\n", "mistakes in calculation. However, there is no need to do this for simpler problems. \n", "\n", "I also **strongly** recommend creating new variables for each problem, for example Rx1, Rx2, etc. for\n", "the range of the random variable in problems 1, 2, etc. That way, you won't have problems if you\n", "forget and use the wrong variable! You can also refer to previous results without problems. \n", "\n", "Following Problem One is an example of what I mean (it is a simple problem, but I am showing you\n", "how you could approach it). \n", "\n", "You are not **required** to do it this way, but I *encourage* you to do something similar. \n", "\n", "[Optional, but a great idea: Write a function `print_solution(...)` which prints out all the answers as shown in the Example Problem!]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example Problem\n", "\n", "\n", "*Describe* the random variable X = \"the number of heads showing on 2 flipped fair coins\"\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Solution:\n", "\n", "(i) Rx = [0, 1, 2]\n", "(ii) fx = [0.25, 0.5, 0.25]\n", "(iii)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAXQ0lEQVR4nO3de7RdZX3u8e9DkIsGQQELhkBQUU7UU2wj1mMrSm2LN9BqLdRbUA719KSodFiwRWqxVqXeasVqVITiBag6asR4GEqFU1Qs4XLEQNGYcgkXTVQKOKpc/J0/5tx0sVh777WTzL2zmd/PGGvseXnXO39r7WQ+e75zrTlTVUiS+mu7uS5AkjS3DAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0DbjCRPSHJFkjuSHDfX9cyWJJXkcdtyHUmWJ7l4tmvS7DAINKUk1yX5zyR3JvlBkk8kWdiuuzDJz9od9+1JLktyYpIdB57/1iR3t8+fePzpJJv7U+DCqtqlqj6wFWofte3btrTfbUW7c763fV23J7kyyQvmui7NPwaBxvHCqloI/ArwVOCkgXUrqmoXYG/gT4AjgdVJMtDmnKpaOPA4dZLt7Aes3ZwCk2w/yarhbe+2Of1vw77Z/m52Az4OnJvkkcONpnh/JINA46uqm4AvA08ase6nVXUhcDjwdOD5M+k7yT8DzwY+2P6F+/gkuyb5hyQbk1yf5KQk27Xtlyf5epL3Jfkx8NYZbu9/JNmUZHE7/8tJbktyYDt/YpLvt0c7Vyd58cBzB7d9W5L1bX/Lk9yY5IdJXj3Q/owkH07ylba/i5LsN0ldOyZ5d5Ib2iOwDyfZebrXU1W/AE4HdgYek+RZSTYkOSHJrcAn2v7/Z5J1SX6cZFWSRw919bz29WxK8jcT7/eIOg9sX8+Pk1yb5GVDr/dDSb7c/i6/nmSvJO9P8pMk/5bkKdO9Js0eg0Bja3eazwOumKxNVd0ArAF+YyZ9V9WhwL/QHGEsrKrvAn8H7Ao8BjgEeBVw9MDTngasBx4FvH2G2/sG8BHgzHZHexZwUlX9W9vk++1r2BX4S+CTSfYe2va3gd2BTwNn0xwtPQ54BU2gLRxo/3LgbcAewJXApyYp7V3A44GD2r4WASdP93rav/iPAe4Evtcu3gt4JM2R1rFJDgXeAbyM5gju+rbuQS8GltEc/R0BvGbEth4GfKV93Y8CjgI+lOSJA81eRnPkuAfwc+CbwOXt/GeB9073mjSLqsqHj0kfwHU0O5fbaHYcHwJ2btddCBwz4jlnAx9tp98K3NU+f+Lx6Em2dV9/wAKaHcjSgfV/SHMOAWA5cMM0tY/a9tcG1j8EuAy4Cvg/QKbo60rgiIFtf29g3ZOBAn5pYNmPgIPa6TOAswfWLQTuBRa380Wz0w/wU+CxA22fDvz7JDUtB+5pX9cm4BLgOe26Z7WvfaeB9h8HTh2q425gyUAdhw2s/yPggoFtXdxO/z7wL0O1fAT4i4HX+9GBdX8MXDP0ft021/+2ffzXw3FDjeNFVfXVGbRfBHxjYP7cqnrFDLe5B7ADTfhMuL7te8KNY/Qz6bar6u4kZwAfAI6vdi8FkORVwPHAknbRwramCT8YmP7Ptr/hZYNHBPfVWlV3tsNZjx56DXsCDwUuGzjFEppQnMwlVfXrk6zbWFU/G5h/NM1f5YN1/IjmPb1uuE6a93t46AiaI4ynDZ14357mqGrC8Hsx1XujOebQkLaqdvjoV2mGebbEJpq/VgfH0vcFbhqY36JL5yZZBPwFzfj5eyY+7dSO338UWAHsXs0J5u/Q7JQ31+KB7S6kGbK5eajNJpqd5BOrarf2sWs1J4M3x/D7czMD72c7xLM7939PFw9M7zuiRmjC4qKBGnerZjjvf21mnZpjBoG2iiQPTXII8AXgX4HVW9JfVd0LnAu8Pcku7c75eOCTW1ws0H6q6Qya4ZLXArfQjOEDPIxmJ7qxbXs0I06Qz9Dzkvx6kh3a7Xyrqu53RFPNCd+PAu9L8qh224uS/M4WbnvCp4GjkxzUht5ft3VcN9DmTUke0Qb664FzRvRzHvD4JK9M8pD28dQk/20r1alZZhBoS30wyR00h/7vBz5HM878i63Q9x/TjJmvBy6m2ZGdPsM+fj/3/x7Bne1O9jjgl4C3tENCR9PsJH+jqq4G3kNzgvMHNGPaX9/C1/JpmqOPH9McMb18knYnAOuAS5LcDnwVeMIWbhuAqroAeAvN7+gW4LE0H/cd9AWa8yZXAl+iCcrhfu4Afrt97s3ArTQnuXccbqv5IQPDopI60J6H2FBVJ03XVpoLHhFIUs8ZBJLUcw4NSVLPeUQgST03775Qtscee9SSJUvmugxJmlcuu+yyTVW156h18y4IlixZwpo1a+a6DEmaV5JcP9k6h4YkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6rlOgyDJYe39TNclOXHE+uVp7kd7Zfs4pst6JEkP1Nn3CJIsAE4DfgvYAFyaZFV7id9B51TViq7qkCRNrcsjgoOBdVW1vqruormP7REdbk+StBm6/GbxIu5//9MNwNNGtHtJkmcC3wXeOHzXJoAkxwLHAuy7774dlKq5tPc++3LrTePcflizZa9Fi7llww1zXYZmSZdBMOr+rsOXOv0i8Jmq+nmS1wFnAoc+4ElVK4GVAMuWLfNyqQ8yt950I/udcN5cl6EB17/rBXNdgmZRl0NDG7j/jbD3YehG2FX1o6r6eTv7UZpb+EmSZlGXQXApcECS/dsbdh8JrBpskGTvgdnDgWs6rEeSNEJnQ0NVdU+SFcD5wALg9Kpam+QUYE1VrQKOS3I4cA/NTb2Xd1WPJGm0Ti9DXVWrgdVDy04emH4z8OYua5AkTc1vFktSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPddpECQ5LMm1SdYlOXGKdi9NUkmWdVmPJOmBOguCJAuA04DnAkuBo5IsHdFuF+A44Ftd1SJJmlyXRwQHA+uqan1V3QWcDRwxot3bgFOBn3VYiyRpEl0GwSLgxoH5De2y+yR5CrC4qs6bqqMkxyZZk2TNxo0bt36lktRjXQZBRiyr+1Ym2wHvA/5kuo6qamVVLauqZXvuuedWLFGS1GUQbAAWD8zvA9w8ML8L8CTgwiTXAb8GrPKEsSTNri6D4FLggCT7J9kBOBJYNbGyqv6jqvaoqiVVtQS4BDi8qtZ0WJMkaUhnQVBV9wArgPOBa4Bzq2ptklOSHN7VdiVJM7N9l51X1Wpg9dCykydp+6wua5EkjeY3iyWp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSem6sIEjyuSTPT2JwSNKDzLg79r8H/gD4XpJ3Jjmww5okSbNorCCoqq9W1cuBXwGuA76S5BtJjk7ykC4LlCR1a+yhniS7A8uBY4ArgL+lCYavdFKZJGlWbD9OoySfBw4EzgJeWFW3tKvOSbKmq+IkSd0b94jgY1W1tKreMRECSXYEqKplkz0pyWFJrk2yLsmJI9a/LslVSa5McnGSpZv1KiRJm23cIPirEcu+OdUTkiwATgOeCywFjhqxo/90VT25qg4CTgXeO2Y9kqStZMqhoSR7AYuAnZM8BUi76uHAQ6fp+2BgXVWtb/s6GzgCuHqiQVXdPtD+YUDNqHpJ0hab7hzB79CcIN6H+/+1fgfwZ9M8dxFw48D8BuBpw42S/G/geGAH4NBp+pQkbWVTBkFVnQmcmeQlVfW5GfadEcse8Bd/VZ0GnJbkD4CTgFc/oKPkWOBYgH333XeGZUiSpjLd0NArquqTwJIkxw+vr6qpxvQ3AIsH5vcBbp6i/dk0X1x7gKpaCawEWLZsmcNHkrQVTXey+GHtz4XALiMeU7kUOCDJ/kl2AI4EVg02SHLAwOzzge+NWbckaSuZbmjoI+3Pv5xpx1V1T5IVwPnAAuD0qlqb5BRgTVWtAlYkeQ5wN/ATRgwLSZK6Nd3Q0AemWl9Vx02zfjWwemjZyQPTrx+jRklSh6b71NBls1KFJGnOjPOpIUnSg9h0Q0Pvr6o3JPkioz/6eXhnlUmSZsV0Q0NntT/f3XUhkqS5Md3Q0GXtz4vaj4AeSHNkcG1V3TUL9UmSOjbuZaifD3wY+D7NN4b3T/KHVfXlLouTJHVvrCAA3gM8u6rWASR5LPAlwCCQpHlu3MtQ/3AiBFrrgR92UI8kaZZN96mh320n1yZZDZxLc47g92guISFJmuemGxp64cD0D4BD2umNwCM6qUiSNKum+9TQ0bNViCRpboz7qaGdgNcCTwR2mlheVa/pqC5J0iwZ92TxWcBeNHcsu4jm3gJ3dFWUJGn2jBsEj6uqtwA/ba8/9Hzgyd2VJUmaLeMGwd3tz9uSPAnYFVjSSUWSpFk17hfKViZ5BPAWmruMLWyn55W999mXW2+6ca7LkLZ9Cx5CMuq245pLey1azC0bbtjq/Y4VBFX1sXbyIuAxW72KWXLrTTey3wnnzXUZGnL9u14w1yVo2L13+39lG9TV/5WxhoaS7J7k75JcnuSyJO9PsnsnFUmSZtW45wjOprmkxEuAlwKbgHO6KkqSNHvGPUfwyKp628D8XyV5URcFSZJm17hHBF9LcmSS7drHy2iuPipJmuemu+jcHTQXmQtwPPDJdtV2wJ3AX3RanSSpc9Nda2iX2SpEkjQ3xj1HQJLDgWe2sxdWlZ8tk6QHgXE/PvpO4PXA1e3j9e0ySdI8N+4RwfOAg6rqFwBJzgSuAE7sqjBJ0uwY91NDALsNTO+6tQuRJM2NcY8I3gFckeRrNJ8geibw5s6qkiTNmmmDIM2Vpy4Gfg14Kk0QnFBVt3ZcmyRpFkwbBFVVSf6pqn6V5sqjkqQHkXHPEVyS5KmdViJJmhPjniN4NvC6JNcBP6UZHqqq+u9dFSZJmh3jBsFzO61CkjRnphwaSrJTkjcAbwIOA26qqusnHtN1nuSwJNcmWZfkAd85SHJ8kquTfDvJBUn22+xXIknaLNOdIzgTWAZcRXNU8J5xO06yADitfd5S4KgkS4eaXQEsa4eYPgucOm7/kqStY7qhoaVV9WSAJB8H/nUGfR8MrKuq9e3zzwaOoLlEBQBV9bWB9pcAr5hB/5KkrWC6I4K7Jyaq6p4Z9r0IGLxT/IZ22WReC3x51IokxyZZk2TNxo0bZ1iGJGkq0x0R/HKS29vpADu38xOfGnr4FM/NiGU1smHyCpohqENGra+qlcBKgGXLlo3sQ5K0eaa7H8GCLeh7A7B4YH4f4ObhRkmeA/w5cEhV/XwLtidJ2gwzuejcTF0KHJBk/yQ7AEcy9M3kJE8BPgIcXlU/7LAWSdIkOguC9pzCCuB84Brg3Kpam+SU9iY3AH8DLAT+McmVSbyEhSTNsrHvULY5qmo1sHpo2ckD08/pcvuSpOl1OTQkSZoHDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSe6zQIkhyW5Nok65KcOGL9M5NcnuSeJC/tshZJ0midBUGSBcBpwHOBpcBRSZYONbsBWA58uqs6JElT277Dvg8G1lXVeoAkZwNHAFdPNKiq69p1v+iwDknSFLocGloE3Dgwv6FdNmNJjk2yJsmajRs3bpXiJEmNLoMgI5bV5nRUVSurallVLdtzzz23sCxJ0qAug2ADsHhgfh/g5g63J0naDF0GwaXAAUn2T7IDcCSwqsPtSZI2Q2dBUFX3ACuA84FrgHOram2SU5IcDpDkqUk2AL8HfCTJ2q7qkSSN1uWnhqiq1cDqoWUnD0xfSjNkJEmaI36zWJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknqu0yBIcliSa5OsS3LiiPU7JjmnXf+tJEu6rEeS9ECdBUGSBcBpwHOBpcBRSZYONXst8JOqehzwPuBdXdUjSRqtyyOCg4F1VbW+qu4CzgaOGGpzBHBmO/1Z4DeTpMOaJElDUlXddJy8FDisqo5p518JPK2qVgy0+U7bZkM7//22zaahvo4Fjm1nnwBc20nR88sewKZpW2k2+TvZNvl7aexXVXuOWrF9hxsd9Zf9cOqM04aqWgms3BpFPVgkWVNVy+a6Dv0XfyfbJn8v0+tyaGgDsHhgfh/g5snaJNke2BX4cYc1SZKGdBkElwIHJNk/yQ7AkcCqoTargFe30y8F/rm6GquSJI3U2dBQVd2TZAVwPrAAOL2q1iY5BVhTVauAjwNnJVlHcyRwZFf1PAg5VLbt8XeybfL3Mo3OThZLkuYHv1ksST1nEEhSzxkE89B0l+7Q7EpyepIftt+L0TYgyeIkX0tyTZK1SV4/1zVtyzxHMM+0l+74LvBbNB+/vRQ4qqquntPCeizJM4E7gX+oqifNdT2CJHsDe1fV5Ul2AS4DXuT/k9E8Iph/xrl0h2ZRVf1f/P7LNqWqbqmqy9vpO4BrgEVzW9W2yyCYfxYBNw7Mb8B/4NKk2qsaPwX41txWsu0yCOafsS7LIQmSLAQ+B7yhqm6f63q2VQbB/DPOpTuk3kvyEJoQ+FRVfX6u69mWGQTzzziX7pB6rb2c/ceBa6rqvXNdz7bOIJhnquoeYOLSHdcA51bV2rmtqt+SfAb4JvCEJBuSvHauaxLPAF4JHJrkyvbxvLkualvlx0clqec8IpCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCNQrSfZJ8oUk30vy/SR/234fY6rn/Nls1SfNBYNAvdF+yejzwD9V1QHA44GFwNuneapBoAc1g0B9cijws6r6BEBV3Qu8EXhNkj9K8sGJhknOS/KsJO8Edm6/kPSpdt2rknw7yf9Lcla7bL8kF7TLL0iyb7v8jCR/314bf32SQ9r7F1yT5IyB7f12km8muTzJP7bXyCHJO5Nc3fb77ll6n9Qznd28XtoGPZHmuvT3qarbk9zAJP8XqurEJCuq6iCAJE8E/hx4RlVtSvLItukHae5HcGaS1wAfAF7UrnsETQgdDnyR5luvxwCXJjmI5vpRJwHPqaqfJjkBOL4NphcDB1ZVJdltK70P0v0YBOqTMPpKrZMtH+VQ4LNVtQmgqibuQ/B04Hfb6bOAUwee88V2R34V8IOqugogyVpgCc2FA5cCX29Gr9iB5pIVtwM/Az6W5EvAeWPWKM2IQaA+WQu8ZHBBkofTXM31P7j/UOlOk/QxbmgMtvl5+/MXA9MT89sD9wJfqaqjHrCx5GDgN2kuLriCJoikrcpzBOqTC4CHJnkV3Hfbz/cAZwDrgYOSbJdkMc2d4Cbc3V7SeKKPlyXZve1jYmjoGzQ7a4CXAxfPoK5LgGckeVzb50OTPL49T7BrVa0G3gAcNKNXK43JIwL1Rjs882LgQ0neQvOH0GqaTwXdBfw7cBXwHeDygaeuBL6d5PKqenmStwMXJbkXuAJYDhwHnJ7kTcBG4OgZ1LUxyXLgM0l2bBefBNwBfCHJTjRHIm/cvFcuTc2rj0pSzzk0JEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HP/HxfGhN3IUL2JAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "Rx0 = list(range(3))\n", "\n", "def f0(k): # if you write f as a Python function, then you can create the list\n", " return C(2,k)/4 # fx by using a list comprehension, as shown here:\n", " \n", "fx0 = [ f0(k) for k in Rx0 ]\n", "\n", "def print_solution(Rx,fx):\n", " print(\"Solution:\\n\")\n", " print(\"(i) Rx =\",Rx)\n", " print(\"(ii) fx =\",round4_list(fx)) # in case you get complicated decimals, round to 4 places\n", " print(\"(iii)\")\n", " draw_distribution( Rx, fx, title='PDF for Example Problem')\n", " \n", "print_solution(Rx0,fx0)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem One\n", "\n", "Do problem 1 from the End of Chapter Problems for Chapter 3 (section 3.3). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution:**\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Two\n", "\n", "Suppose you deal a 5-card hand from a standard deck which has been shuffled well. \n", "\n", "Let $X$ = \"The number of Spades occurring in the hand.\" \n", "\n", "*Describe* the random variable $X$. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# template for solution\n", "\n", "# Range\n", "Rx2 = [] # your code here\n", "\n", "# define f(k) = P(X = k)\n", "\n", "def f2(k):\n", " pass # your code here\n", "\n", "# PMF\n", "fx2 = [ ] # your code here\n", "\n", "# if you want to write print_solution\n", "\n", "#print_solution(Rx2,fx2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Three\n", "\n", "This random experiment takes the form of a game, in which one \"trial\" consists of drawing two chips at random and without replacement from an urn that contains 2 red, 4 white, and 5 blue chips. For\n", "each blue chip we win $\\$$1, for each white chip we win $\\$$2, but for each red chip we \n", "lose $\\$$3. \n", "\n", "Describe the random variable $Y$ which implements this random experiment (a \"win\" returns a positive number, whereas a \"loss\" returns a negative number, e.g., -3). \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution:**\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Four\n", "\n", "We refer to the random variable $X$ from Problem Two. \n", "\n", "Describe the random variable $Y = 2X + X - 1$\n", "\n", "\n", "Hint: when more than one instance of a random variable is involved, it is often useful to draw a matrix of all possibilities. Consider the two instances\n", "of $X$ (i.e., two \"pokes\" at the same random variable) and draw a matrix of each of the possible outputs of the random variables, one along the rows and one along the columns.\n", "\n", "| | 0 | 1 | 2 | 3 | 4 | 5 |\n", "|----|-----|-----|-----|-----|-----|-----|\n", "| 0 | | | | | | | \n", "| 1 | | | | | | | \n", "| 2 | | | | | | |\n", "| 3 | | | | | | | \n", "| 4 | | | | | | | \n", "| 5 | | | | | | | \n", "\n", "For each slot in the resulting matrix, you can figure out the probabilities\n", "by multiplication (since the two instances of X are independent), as shown in class. The value in each slot will be $2X + X - 1$. \n", "Some outcomes will be the same, and you will have to add the probabilities\n", "to get the PMF of Y. \n", "\n", "I would do this calculation in Python, but it is up to you. \n", "If you use Python, you might want to use a dictionary to keep track of the values in the PMF. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Five (Special Distributions)\n", "\n", "\n", "Suppose numbers in the range $[0..1)$ are randomly and independently selected with replacement\n", "and rounded to 3 decimal places. Therefore we can assume that all possible combinations of 3 digits after the decimal point are equally likely. \n", "\n", "For each of these, specify precisely the discrete distribution involved, including all parameters; phrase it as a probability in the form `P(...X....)` and answer the question. \n", "\n", "(a) What is the probability that the first selection is no more than 0.345?\n", "\n", "(b) What is the probability that 0.345 occurs at least twice in the first 1000 selections?\n", "\n", "(c) What is the probability that 0.345 is selected for the first time on the 1000th selection?\n", "\n", "(d) What is the probability that 0.345 is selected for the fifth time precisely on the 1000th selection?\n", "\n", "Hint: This uses all 4 of the discrete distributions we studied in Lecture 11. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution:** \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Six (Binomial)\n", "\n", "Answer the following questions, each of which involves an appropriate binomial distribution. \n", "\n", "(A) If you roll a die 5 times, what is the probability that on exactly two of these rolls you get either a 5 or a 6?\n", "\n", "(B) If two fair dice are rolled 10 times, what is the probability of at least one 6 (on either die) in exactly five of these 10 rolls?\n", "\n", "(C) Suppose that each day the price of a stock moves up 1/8th of a point with probability 1/3 and moves down 1/8th of the point with probability 2/3. If the price fluctuations from day to day are independent and identically distributed, what is the probability that after 6 days the stock has its original price?\n", "\n", "Hint for (C): you could draw a tree, but it is easier to answer the following question: how many moves up and how many moves down result in no change in the stock price after 3 days?\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution:** \n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Seven (Binomial)\n", "\n", " A restaurant serves 8 fish entrees, 12 beef, and 10 chicken. If customers select from these entrees randomly and independently, and the kitchen has plenty of each entree (so this is \"with replacement\"), what is the probability\n", "\n", "(A) that at least 4 of the next 10 customers order beef?\n", "\n", "(B) that between 2 and 6 (inclusive) of the next 20 customers order fish?\n", "\n", "(C) that of the next 10 customers, at most 1 of the first 5 orders chicken, and at least 2 of the last 5 orders beef. \n", "\n", "(D) that of the next 5 customers, the first orders chicken, the second orders fish, and of the remaining 3, at least 2 order beef?\n", "\n", "**NOTE: For each of these, you must state the exact distribution being used, with all parameters, and give the solution in the form:** \n", "\n", " X ~ (exact distribution)\n", " P(...X...) = (formula) = (numeric answer). \n", "\n", "Hint: For (C) and (D) remember that all customers choose their meal independently. \n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution:**\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Eight (Binomial)\n", "\n", "Suppose a professor of probability is tired of reading the depressing news and so he decides that he will quickly scan the first 5 headlines in the New Yorks Times and the first 5 headlines in the Boston Globe and if at most 3 of the articles in each are depressing, he will read the news that day. Further suppose that the probability of a NYTs headline being depressing is 0.6 and for the Globe the probability of a headline being depressing is 0.55. \n", "\n", "(a) What is the probability that he will read the news the first day he tries this?\n", "\n", "(b) In order to be \"well-informed\" he needs to read the news at least half the time; what is the probability that he will be well-informed after doing this for a week? \n", "\n", "Hint: This is another problem where there are two independent parts of the random experiment.\n", "You might want to phrase it as three different random variables, all three being binomial.\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution:**\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lab Instructions\n", "\n", "In previous homeworks, we investigated how to generate random values (called \"random variates\" in the literature). In this lab, we will make this more explicit by considering three different ways for simulating a random variable:\n", "\n", " > 1. By simulating the original physical experiment;\n", " > 2. By inverting the CDF; and\n", " > 3. By using an explicit formula.\n", "\n", "In the last two, we will convert a random variate in the range $[0..1)$ created by\n", "the `numpy` function `random()` into a random variate from a different distribution. \n", "You may not use any other function from the `numpy` random library for the following problems. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Nine: Generating a Binomial Distribution by Simulation\n", "\n", "For this problem, complete the function stub to simulate the\n", "a random variable $Y\\sim B(N,p)$, which describes a situation in which\n", "we have a (possibly) unfair coin where\n", "the probability of heads is $p$ and we flip this coin $N$ times. \n", "\n", "We create random variables $X$ and $Y$ to solve this, as shown in lecture. \n", "Let $X\\sim \\text{bernoulli}(p)$ and $Y = X_1+X_2+\\cdots + X_N$. \n", "\n", "Demonstrate your solution by generating $10^5$ random variates for the parameters shown and displaying them using show_distribution(...).\n", "\n", "Note: You MUST use `random()` to implement this, NOT `randint(2)`. Think about using a list comprehension and `sum` to implement the binomial. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# First, create a function which simulates the coin, where\n", "# you return 1 with probability p and 0 with probability 1-p. \n", "\n", "num_trials = 10**5\n", "\n", "N9 = 10\n", "p9 = 0.25\n", "\n", "seed(0)\n", "\n", "# Simulate the random variable X ~ bernoulli(0.25)\n", "\n", "def bernoulli(p):\n", " pass # your code here, you already wrote this in previous hws!\n", "\n", "def X(): \n", " pass # your code here, just call the previous with the appropriate parameter\n", " \n", "# simulate the random variable Y ~ B(10,0.25)\n", "\n", "def B(N,p):\n", " pass # your code here\n", "\n", "def Y():\n", " pass # your code here, just call the previous with the appropriate parameter\n", "\n", "\n", "# Now show the result of 10^5 trials with show_distribution\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Ten: Generating a Distribution by Inverting the CDF\n", "\n", "In this problem we will investigate how to implement a random variable given by an arbitary probability distribution function. First, however, you need to know about the CDF, the Cumulative Distribution Function, which is a simple but crucially important idea, which you can understand by looking at the last slide on Lecture 11 or here. Please look at this before continuing with the lab. \n", "\n", "We will return to an idea (sometimes called the \"pinwheel method\") we used in a previous homework: we can essentially \"invert\" the CDF to obtain a function from a random variate in the range $[0..1)$ into a random variable from the given distribution. \n" ] }, { "attachments": { "Screen%20Shot%202021-03-12%20at%2011.00.35%20AM.png": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAFCCAYAAACQFIg4AAAK12lDQ1BJQ0MgUHJvZmlsZQAASImVlwdUU9kWhs+9N52EFoiAlNA70gkgJfTQexOVkAQSSgwJAUTsDI7gWFARAXVER0UUHB2KjAWxYBsUe58gg4I6DhZsqMwNPMLMvPXeW2+vdXK+tbPPPnufe27WHwAoIWyRKBdWBSBPWCCODfajJ6ek0vGDAAEYQAGeAMfmSETM6OhwgNrU/Hd7dwtA8vm6rTzXv3//X02dy5NwAIDSUM7gSjh5KHeh4xlHJC4AADmA+o2LCkRyvoqyhhgtEOXf5Jw1yR/knDHBGPJETHysP8p0AAhkNlucBQDZBvXTCzlZaB6yvAd7IVcgRLkUZW8On81F+RjKNnl5C+Q8hLIFGi8CgIKeDmBk/CVn1t/yZyjys9lZCp7sa8IIAQKJKJe98P88mv9tebnSqT3M0EHmi0Ni5fuh53cnZ0GYgoUZkVFTLOBO1iRnvjQkYYo5Ev/UKeayA8IUa3Mjw6c4UxDEUuQpYMVPMU8SGDfF4gWxir0yxf7MKWaLJ/YloSyT5iQo/HweS5G/hB+fNMWFgsTIKZbkxIVNx/gr/GJprKJ+njDYb3rfIEXveZK/9CtgKdYW8ONDFL2zp+vnCZnTOSXJitq4vIDA6ZgERbyowE+xlyg3WhHPyw1W+CWFcYq1BejlnF4brTjDbHZo9BQDAYgAbMChq0wRAAW84gJ5I/4LRAvFgix+AZ2Jvm08OkvIsbOhO9o7OgAgf3cnr8Mb2sQ7CdEuTvtWFaFXvRUF+rQvRB+AIwHoYxmZ9pn7AqCMB+C8HkcqLpz0YeQfWPTpqQANoA30gTGwALbAEbiivxG+IBCEgigQD1LAPLRWPsgDYlAESsFyUA4qwXqwGdSCHWAX2AcOgsOgHRwDp8A5cAlcBTfBfSADg+A5GAHvwBgEQXiIAlEhbcgAMoWsIUeIAXlDgVA4FAulQOlQFiSEpFAptBKqhKqgWmgn1Aj9CB2FTkEXoD7oLtQPDUOvoU8wApNhDVgPNoNnwQyYCYfB8fBcOAvOh0vgMngtXAM3wAfgNvgUfAm+Ccvg5/AoAhAlhIYYIrYIA/FHopBUJBMRI0uQCqQaaUCakU6kB7mOyJAXyEcMDkPF0DG2GE9MCCYBw8HkY5Zg1mBqMfswbZgzmOuYfswI5iuWgtXFWmM9sCxsMjYLW4Qtx1Zj92BbsWexN7GD2Hc4HI6GM8e54UJwKbhs3CLcGtw2XAuuC9eHG8CN4vF4bbw13gsfhWfjC/Dl+K34A/iT+Gv4QfwHghLBgOBICCKkEoSEFYRqwn7CCcI1wlPCGFGVaEr0IEYRucSFxHXE3cRO4hXiIHGMpEYyJ3mR4knZpOWkGlIz6SzpAemNkpKSkZK7UoySQGmZUo3SIaXzSv1KH8nqZCuyPzmNLCWvJe8ld5Hvkt9QKBQzii8llVJAWUtppJymPKJ8UKYq2ymzlLnKS5XrlNuUrym/VCGqmKowVeaplKhUqxxRuaLyQpWoaqbqr8pWXaJap3pU9bbqqBpVzUEtSi1PbY3afrULakPqeHUz9UB1rnqZ+i710+oDVIRqTPWncqgrqbupZ6mDGjgNcw2WRrZGpcZBjV6NEU11TWfNRM1izTrN45oyGkIzo7FoubR1tMO0W7RPM/RmMGfwZqye0Tzj2oz3WjO1fLV4WhVaLVo3tT5p07UDtXO0N2i3az/UwehY6cToFOls1zmr82KmxkzPmZyZFTMPz7ynC+ta6cbqLtLdpXtZd1RPXy9YT6S3Ve+03gt9mr6vfrb+Jv0T+sMGVANvA4HBJoOTBs/omnQmPZdeQz9DHzHUNQwxlBruNOw1HDMyN0owWmHUYvTQmGTMMM403mTcbTxiYmASYVJq0mRyz5RoyjDlm24x7TF9b2ZulmS2yqzdbMhcy5xlXmLeZP7AgmLhY5Fv0WBxwxJnybDMsdxmedUKtnKx4lvVWV2xhq1drQXW26z7bLA27jZCmwab27ZkW6ZtoW2Tbb8dzS7cboVdu93LWSazUmdtmNUz66u9i32u/W77+w7qDqEOKxw6HV47WjlyHOscbzhRnIKcljp1OL1ytnbmOW93vuNCdYlwWeXS7fLF1c1V7NrsOuxm4pbuVu92m6HBiGasYZx3x7r7uS91P+b+0cPVo8DjsMcfnraeOZ77PYdmm8/mzd49e8DLyIvttdNL5k33Tvf+3lvmY+jD9mnweexr7Mv13eP7lGnJzGYeYL70s/cT+7X6vff38F/s3xWABAQHVAT0BqoHJgTWBj4KMgrKCmoKGgl2CV4U3BWCDQkL2RBym6XH4rAaWSOhbqGLQ8+EkcPiwmrDHodbhYvDOyPgiNCIjREPIk0jhZHtUSCKFbUx6mG0eXR+9M8xuJjomLqYJ7EOsaWxPXHUuPlx++PexfvFr4u/n2CRIE3oTlRJTEtsTHyfFJBUlSRLnpW8OPlSik6KIKUjFZ+amLondXRO4JzNcwbTXNLK027NNZ9bPPfCPJ15ufOOz1eZz55/JB2bnpS+P/0zO4rdwB7NYGXUZ4xw/DlbOM+5vtxN3GGeF6+K9zTTK7MqcyjLK2tj1jDfh1/NfyHwF9QKXmWHZO/Ifp8TlbM3Zzw3Kbclj5CXnndUqC7MEZ5ZoL+geEGfyFpULpLle+Rvzh8Rh4n3SCDJXElHgQYqki5LLaTfSPsLvQvrCj8UJRYdKVYrFhZfXmi1cPXCpyVBJT8swiziLOouNSxdXtq/mLl45xJoScaS7qXGS8uWDi4LXrZvOWl5zvJfVtivqFrxdmXSys4yvbJlZQPfBH/TVK5cLi6/vcpz1Y5vMd8Kvu1d7bR66+qvFdyKi5X2ldWVn9dw1lz8zuG7mu/G12au7V3num77etx64fpbG3w27KtSqyqpGtgYsbFtE31Txaa3m+dvvlDtXL1jC2mLdIusJrymY6vJ1vVbP9fya2/W+dW11OvWr65/v4277dp23+3NO/R2VO749L3g+zs7g3e2NZg1VO/C7Src9WR34u6eHxg/NO7R2VO558te4V7Zvth9ZxrdGhv36+5f1wQ3SZuGD6QduHow4GBHs23zzhZaS+UhcEh66NmP6T/eOhx2uPsI40jzT6Y/1bdSWyvaoLaFbSPt/HZZR0pH39HQo92dnp2tP9v9vPeY4bG645rH150gnSg7MX6y5ORol6jrxamsUwPd87vvn04+feNMzJnes2Fnz58LOne6h9lz8rzX+WMXPC4cvci42H7J9VLbZZfLrb+4/NLa69rbdsXtSsdV96udfbP7TlzzuXbqesD1czdYNy7djLzZdyvh1p3babdld7h3hu7m3n11r/De2P1lD7APKh6qPqx+pPuo4VfLX1tkrrLj/QH9lx/HPb4/wBl4/pvkt8+DZU8oT6qfGjxtHHIcOjYcNHz12Zxng89Fz8delP+u9nv9S4uXP/3h+8flkeSRwVfiV+Ov17zRfrP3rfPb7tHo0Ufv8t6Nva/4oP1h30fGx55PSZ+ejhV9xn+u+WL5pfNr2NcH43nj4yK2mD0hBRB0wJmZALzei2rjFACoqC4nzZnU1hMGTf4fmCDwn3hSf0+YKwDNqOaIQZGJzkfkchadKeiQS6J4XwA7OSnGv0yS6eQ4mYuMKkvsh/HxN3oA4DsB+CIeHx/bNj7+ZTda7F0AuvInNb3ccKiWbzbQGxiOvlH8K/inTer9v/T4zxnIK3AG/5z/BO/lFXZGZwNkAAAAlmVYSWZNTQAqAAAACAAFARIAAwAAAAEAAQAAARoABQAAAAEAAABKARsABQAAAAEAAABSASgAAwAAAAEAAgAAh2kABAAAAAEAAABaAAAAAAAAAJAAAAABAAAAkAAAAAEAA5KGAAcAAAASAAAAhKACAAQAAAABAAABsKADAAQAAAABAAABQgAAAABBU0NJSQAAAFNjcmVlbnNob3RVH5/fAAAACXBIWXMAABYlAAAWJQFJUiTwAAACc2lUWHRYTUw6Y29tLmFkb2JlLnhtcAAAAAAAPHg6eG1wbWV0YSB4bWxuczp4PSJhZG9iZTpuczptZXRhLyIgeDp4bXB0az0iWE1QIENvcmUgNi4wLjAiPgogICA8cmRmOlJERiB4bWxuczpyZGY9Imh0dHA6Ly93d3cudzMub3JnLzE5OTkvMDIvMjItcmRmLXN5bnRheC1ucyMiPgogICAgICA8cmRmOkRlc2NyaXB0aW9uIHJkZjphYm91dD0iIgogICAgICAgICAgICB4bWxuczpleGlmPSJodHRwOi8vbnMuYWRvYmUuY29tL2V4aWYvMS4wLyIKICAgICAgICAgICAgeG1sbnM6dGlmZj0iaHR0cDovL25zLmFkb2JlLmNvbS90aWZmLzEuMC8iPgogICAgICAgICA8ZXhpZjpVc2VyQ29tbWVudD5TY3JlZW5zaG90PC9leGlmOlVzZXJDb21tZW50PgogICAgICAgICA8ZXhpZjpQaXhlbFlEaW1lbnNpb24+NjI4PC9leGlmOlBpeGVsWURpbWVuc2lvbj4KICAgICAgICAgPGV4aWY6UGl4ZWxYRGltZW5zaW9uPjg0MjwvZXhpZjpQaXhlbFhEaW1lbnNpb24+CiAgICAgICAgIDx0aWZmOk9yaWVudGF0aW9uPjE8L3RpZmY6T3JpZW50YXRpb24+CiAgICAgICAgIDx0aWZmOlJlc29sdXRpb25Vbml0PjI8L3RpZmY6UmVzb2x1dGlvblVuaXQ+CiAgICAgIDwvcmRmOkRlc2NyaXB0aW9uPgogICA8L3JkZjpSREY+CjwveDp4bXBtZXRhPgpoB9lXAABAAElEQVR4Ae2dB5hV1bm/vyn0JqiAgAIiYAFEIiq22EWTGFs0liReY5KbBH2SPEme5BrN/6aZasxNbhJNNMYWYy+5FhBFxY4FCxaMUpSigHSkzfz3+w3fcc/hDHMGpuwz81twZpfV33XO+u1vrbX3LqtOnMmJgAiIgAiIQIkRKC+x8qq4IiACIiACIuAEJGD6IoiACIiACJQkAQlYSTabCi0CIiACIiAB03dABERABESgJAlUbk2p61v3UVZWtjXJNnmcdLkLlTHtT2EKhdlSIfPjF0qjmDBbyqNYv83zKUvqU2zsmnCbp1FzPrjU59+w3LYcuqqqKtcekf+WY2zuu1l5Ex4JlVoB02G2Np9aCepABESgyQiUJT/YRl+FSJJZ//EXU8ZiwjRZy2xDwnWVu9qSdsnrsLchm5KKWieTEviulhRoFVYEmpHAVllgy5Yts+XLl7tIpYWKq+RevXpZly5djH0c/nzoQPiXO1dkR0q8Qi6dbyH//HPr1q2z999/3zZs2GA9e/a07t271wpCPrfffrvddttttscee9iXvvQl6927t5e7vryicyT9f/zjH/bCCy9Y165dPY2RI0fm0lizZo1df/31Nn/+fD/34Ycf2vDhw+3zn/98owl+lIU2+stf/mLPPfecfepTn7LPfOYzVllZmStLrcqnDiL+ihUr7Nprr7X33nvPy1ZeXu5xTzjhBBs9erTHeOedd+yaa66xjRs3+jHfCfwPOeSQevNJZVnv7urVq+21116zhx56yNq3b2/nn39+vXHSAaJOS5cutauuusqmTZtmxx9/vH32s591JoSNMKtWrbJFixZZRUWF7bDDDtaxY8d0UtoXARHIEIEGzYHxI8fRCe+yyy7+GTBggMWHc/fcc4+HocPjE50/2/Ky5FzyaYgVQLxCH8+kiD8hpHRalHPQoEH2+OOPe0z8wv+OO+6wU045xTu2F1980TtJBIa8o96FsouOj052xIgR9ve//907RsRv1KhR9u9//zvHgM4RYaSD7NChg3eSCOqW0i+U55bOUV7c9773PRfjE0880c4880z75z//uaVoOb8oy+zZs+3rX/+6de7c2T/dunWzH/7wh/bWW2/lws6ZM8cuvPBCrwfiyIf6NJaLsrzxxhv27W9/2yZNmmR33313g3nBBJH90Y9+ZLfccovB5HOf+5y3FWUln8iLth+UfEd23nlnmzJlilclviN+oD8iIALZIZD8cBvsEkui+ve//331UUcdVf3uu+9WL1iwoDq5+q7+5S9/6R8STDrA6lmzZlW//fbb1Yn1Uz1v3jw/Tjr06pUrVxaVZ9LpVJNXIiSbfYpKIAlEGriJEyeivl5WyoNLOibfJuJTvddee1XfcMMNfpxYSB722Wef9eNIww/q+EMZE2ukevHixR5i/fr11R//+Mer//CHP+Ri4EcZEgHLnWMnylHr5FYcRDlfeeUVzwf2uAcffNCPaSNcMfnBnXYNl1h01bvuumv13Llz41T11KlTq7/85S/njmOnmPQjbDFbWOKox/jx44uJkgsTZZk5c6YzYIuL78OSJUv8OMLx3aCdEgu6OhF99wuufqA/IiACmSHQIAsM2U1+6D6s0rdvXx/O6devny1cuNASgbKDDjrIr1yTH7x99atf9StZhmqwxO68804/HjJkiHHljkso+LbQH9L4r//6Lxs5YqQde+yxduSRR/qH9BjWeeSRRzwa5SnGUQbcjjvuaO3atauVd9JRW9Lp28EHH5wLc9555/lQEyfCqnHPOv5gUXFVzxAqjiG4hx9+2HbaaadcjKrqmrJiwf7ust8ZVl8i5vVaebkEitwJvv379/cYMMcxHFisg3GfPn1yw4MM3yUib5Em6TDM9vTTT9vVV19tl156qT355JP+/ajPai22DBGOfHBYxFtr4b388suWXFDkyj9mzBhPk2FQHN9FPnw3aEMssPjOeAD9EQERyByBBgtY1IAfO8OFidVle++9t3dkBx54oJ100knesTFU89vf/taH0JKrdkuuon0uhuEz5piIvyVhKK8ot6997Wt297/utiuuuMKuvPJK//zpT3+yV1991YfnKEuxnQz54ShHuDgXnSIihKPDHDx4sFHufEeciJfvh+iGo5y4xEqNU7mh0+kvTLdly5c5K4Se4cf8Tj/yqSuvXKIFdujo6ayDL8OAOObecOk0Ix/3yPuDX4jH5Zdfbuecc46nGbzo7JnvY4iPtMeNG+dzYnnJ+OGW8ikUPn0uyhv1Sfvl79eVz7z581yAO3Xq5FEoOy7qwn6kz0XR2rVrOSUnAiKQYQJbtYiD+tBRbLfddt4JJ0MuubkkRAA/OopvfOMbvqADMTj11FN9UUCPHj3cPzqLOtkkepMM91gyROlXxRGOeHQuhx56qOdPXvWmFZGTbbFh6cSi885FT8pUV3zKEeHp7JkbwqpjsQhpIbTUnbklrBjC/sd//IfPI5577rl2+OGH1+JSVz65smxhh7xYtBKOsuGYo8p3W8qHWMyoMY937733GvXCxUUDC1CwYMIqO+KII3wO8eijj/Zz6bbZUj6eaAP+pNPNj1ZXPhXlFTlrMj9O/nFdaeSH07EIiEDLEti8RyuyPBurNvqQ4dChQ+073/lOzrKJzjKSiU4zmWvwlYvFChjpzJgxw1fRsaIvhgrp+EmLBRMI47a46KjC8mLYj8UXdP5YFkz213JJb85KQhZjMMwUHXm6Q8XywnJkMcCee+7p0SMcLFjogqM+DFPR6cfQXpQHf4YWyYuhPBZQNMQR54knnsgxo144hntx6XxYmYfjYiTfIV64Rx991D75yU/6IhiOqQ91pl34YMVQNywwHBYlLriwJR/CET6sIA/UgD9R7tgWisrFFO1JPmnHwoxf//rXPrQLT9jiou3TYbUvAiJQIgSSzqVBLia0meA++eSTay0ISDrlWsfJsF/1AQcc4As4WCCRIKmOhQWEbQ4X5U1WsHn+ifXm2abzTzrW6mOOOaY6GfJ0PxaaUNZEQP040kiEqzpZReh+yWq4zdK5+m9Xu1+yytH9iJ8Ime/zh/gshgg3K1nkQj7JXJmfinw4+MEPfuB+LAzBpf38RIE/UadkpaDHTUTYQyXDuX68Zs1qP460kqHY6uSCoDoRpOo333zT/SKNpCX9mEUNyXxndXJ7Qa24HFAX6hQuGVL2fFiwg4t8WBQxcOBA95s+fXotPz8o8s99991XnSzTr/UdI2qUOVk16HkQJhZnRBmSeUH3CybxfWQBTzqN2GfBys033+x+kYYf6I8IiEBmCHCVXLSLjiJZkOGdAZ1vci+Nd/z8yMOfVXzJ0m0Pkyy68PSTIahcnOjgI3xdBcB/S5+64qXPR+cTAsZqQVzkHVvKRH1YScn2+9//fq4DjjB0xPjxueSSSzydSJ+OMfySxSzVu+++ux9fdNFFHo4/ISQ//vGPq//85z+7P4IYIhD5sEpxv/32c/9kiM7jh18usXp2kkUVHv+yyy7zLfXHkU6UOVkI436U+4knnnT/8IttcGFVKS7Osz/h6xM8fmJ1VidL1H0/Vu6l82G1Inl885vfrBWfNLbkos6JNenfp2B6xhlnVP/85z+vjtWJEQ6xjzbg4gAXfuxzgYJ/sEEQcekwHFNHCRgk5EQg2wQa9CSOpCo+/MRNpQzVMFTFcBHDNdzcGkM73ETL/A/DMwwb4ceCCOKQBivzkivy3BBT0qk0mUs6Jx/ymjx5si+ooDwMQ0ZdyDj2KTNzPUlH6WGpX/hFAak783KsYuOG6EifITIWMzC8Rh6wIG7cI0d8hgVZDffMM884CxZaJBaqD6kRlg/xCcMN0AxH/ud//udmZYiyFNpGeVmswkpN7n9jTipW3REnHQZ/8sQ/VmdG2dmyYIY6szgn37Fwg7IyXMnQMPN43PsWLtjceOONloiOl+VjH/tYjlmEq2sb5eS781Zy/xlDj5SVRSrc0Lzvvvv6cYSDL/VheBd+8X0M/w3rN9hDUx6y559/3ts3zYQyRDj2Gf7lvjHmbqMenJcTARHIDoEGCRjFTv/IC1WjPv+IU2y4CL+12+h8kvt+fDk+IsocSKxCi3QLlafQuQjPNvxjm/ZryH7Ejy1PwODpHCyHZ54szhebZkPD15VuIqm5lZN1hanrfJQBsUmsSRs2bJjfTM3FQ0NcpFNfnPrC1ecf6SP8lJnFNwjv6aefLgELONqKQMYINHgZfVydIwzpT9SrkD9+dCARvtjOJNJsjC0dEhbT/vvvb1hjOMoTjnJznxZli09cwUcYtlGPdB0K1Tld14ifjss5wqTTiXCs7rv11ltdvDhXqBwRttA2V55N9UnXMz98oXJGGJ6YQvn4FHJRH/zS+xxHmbGUWJF58cUX5yxf/It1ubrkfd/y65QfLj/9nH8dTCI9rON99tnHOnfqnLunLz8tHYuACGSDQIMtsGwUu+GloIPiyjqZ3/Chz3wLrOEpKkZrJMDwLysUGaqMIcvWWE/VSQRaA4E2I2Cl1Fhh8YQVU0plL1RW6tNa6lKofjonAiLQMgTalIC1NmFoma9M689V35PW38aqYesg0KYErHU0mWohAiIgAiIAgQYv4hA2ERABERABEcgCAQlYFlpBZRABERABEWgwAQlYg5EpggiIgAiIQBYISMCy0AoqgwiIgAiIQIMJSMAajEwRREAEREAEskBAApaFVlAZREAEREAEGkxAAtZgZIogAiIgAiKQBQISsCy0gsogAiIgAiLQYAISsAYjUwQREAEREIEsEJCAZaEVVAYREAEREIEGE5CANRiZIoiACIiACGSBgAQsC62gMoiACIiACDSYgASswcgUQQREQAREIAsEJGBZaAWVQQREQAREoMEEJGANRqYITUGAd3DFe7i2Jv364qfTzg+b9tuavLcmTn4ZtiaNdJx0HRo7bfJJp5/OtzH2o7yxbYw0lUbbIKD3gbWNdm6SWhbq1Er9zctRp62tR8TPB74t6W1t3PwyFHscdWjufIstX33hCpW/0Ln60pF/9glIwLLfRq2+hHQuK1assLVr11qPHj2sffv2DarzunXrbNnyZda+XXuPnx95zZo1tnz5cuvWrZt16NDB89qwYYNtv/32RtxFixb5fseOHfOjNuox9UQUKM/KlSutU6dO1rVr123OI+qw3XbbWefOnW3JkiWeZvfu3a2ysnKb0//www9t2bJlXlbSb2xho22WLl1q5eXl1rNnT+vSpcs2l1kJtA0CFf8vcW2jqqplYxJYtWqVvffee97x0AHRIeMQCFxc8fpB8ic6vULnEa67777bbr31Vhs+fLh3YlVVVbk4kUZ+XM6T7ttvv21//N8/ujCNGDGiVt74v/nmm/a3v/3NdthhBxfHW265xR544AE76KCDbMGCBXbRRRfZqFGjDAGgo6Y87dq1s7LysqQiNblH+TmKcqTP1YQy5wEXmASX9evXe3oVFRUejPLceOONXpadd97Z6xBpRjqkjTAtXrzY/SMu59NhOZ4/f779+Mc/tiFDhtiOO+5ol19+uU2bNs2PEW1Y4qK8Eb/QsV8MJAyizIR56623nB+C2K9fPxeaSCNd3tjP90vnHWFiy4UL7XHhhRfaE088YUOHDrX+/ftvVscIX1/a1JU0aUcuSBBF4nDRgLDTro0h6lEebVuWwLZfnrVs+ZV7MxOgM6BTe+qpp+yLX/yi9enTx0vAlfMJJ5xgp556qnei0TnmF6/QeTodOpx33nkn19nS8eS7QnEJg+AQd7fddvMo+eEGDx5sEyZMMDpg8kEUsLqIR2f5s5/9zHr37u2C89e//tUttbPOOsu3lmhYvstPP/xhc9ttt9l3v/tdGzpsqG3csNFFE3E87rjjbOzYsW4hUZ7zzjvPO9hIK7aRFtt58+bZD3/4Q/vyl7/sYht++WFpA65DEWCEB2tm9erVORHIZ5kfP308e/Zsu/baa23vvfe2U045xbPcZZddnB/WV6SVjhPliu2W/CJMbBHzyZMn21e/+lU7+uijc9ZXXWnUdT6+l7Qp6d155532ve99z3bffXcXr3vuucf+9a9/2Q9+8AP/nkT4KIe2pUlAAlaa7dbipcbi4mof64Ur/yeffNL+8pe/eIf9iU98wq0HOhuu6Oks6FyxIrBKOEdHSIfIcBHh8OOKmc6XK2U6YoYTY1iPjok8OU94ht+IS7w4dotlyWJbv269x8P6wD/S5sqbsGxjmBK/6Jgp2xtvvOFCN3fuXLc2CI8jL+JRF67mKQfxuKJPd4aI47HHHusdMhbfu+++a/fdd59z+s53vmOwIR0sVeLiKDd5M6yJo87U/f3337fHH3/cjjnmGBs4cKAzJA4WBuUnHuKPMFMW0t24caNv2f/ggw88H8pHGPLEn/ikE0N1WNOkBS/Yv/TSS75P2TkHa8pEmvAgPdoCDuzj361rN7dY161P6rJsuZdvfbK/IRHxGCoN8fNKJn8QWQSMNqd8UUbYUka2xGGYlTRw5Mt3gfKTP3WCVZSLcHvssYfdftvtbtkhjFiofDcRZNpErvUQ0BBi62nLZqsJncXMmTO98zn99NN92I/O7ZFHHrGddtrJO5w//elP9sILL9gNN9zgV8R77rmnD9dx/n//93/toYce8o5o0KBB3hkR9tVXX7V///vfHochRTojLCQ6TwTyf/7nfzwuV9cIBXkxj8WQ3bPPPutDiU8//bSHo/PFqsI6eeWVV+zXv/619erVyzvJF1980TtqrCKE6oILLjCsohkzZthVV13l5ZgyZYoPndH5MdyIP1YmnSr5I0qICudwIXSUE/epT33KqBvWC3EZ5kQcDzjgAN//wx/+4IJA/Z577jn79a9+bVh/V199tceno7355pvtjjvu8HhYFXTwdNo//elP3eLE/7HHHnOR+c1vfmMMRzLER9mxpOBJfR588CEvJ3khir///e9dCGgT3P333+95wQfL+re//a2HoyzwpW6wpz2oD0yuueYaow633367t+uAAQNcYN94/Q37/ve/bwsXLkyGhG+zP//5z84aBpQ/xJ4trPg+0HaTJk0y5tqoA6JGGeHx8MMPJ9+PSiN9ygGPP/7xj26dXnrppS7gI0eO9IsU6kIY2gRRY5gWoeN7yXfonHPOsb59+xIs115+oD8lS2DzcZqSrYoK3twEGI6jo+JDp0PnSAeBJUEHTwd12mmn+fwGHRadFQJw8cUX22c+8xnvrKdOnWpYAHQ8zIUgOt/4xjd8yO2yyy5za4B60RExREnnxRAdlgmdHlfjxMWCmTVrlh1yyCH23//93y5q//jHPzxtOn06ScKkO1DiEZ982acjPPjggz3vX/ziF3b++ef7UBp1mTNnjuPFakBM6MzpJAs56o9Fg0PYmdcbN26cCxXCShpYOdQbAWVoq1//fkaHDKN99tnHheioo45ya4LhTwTq4x//uJcf8UCUuXigjAjDxIkTXdzIk7pcd911bhnXzO+NdCF5+eWXvW1ef/11t84Ii6MdKQ+dPEOHiO/JJ5/sFxLHH3+8x2FOjfJy4fB///d/LpJYN9/85jeN9MgPi4/2hTsXA7Qx/v/85z+97jChbNEGzFd++tOftr322svb7Nxzz/X0EU4sP+b1jjzyyESE/2aPPvqoW4+U9e9//7uXG3+sUxjjIm2sahhiCf/qV7/y9mJIGBGUa10ENITYutqz2WpDZ4EVxRBWDE0xREMni0XEVfBhhx3mx3SwXNnfe++9flWMFcIcFFYXFgTzFAyJMYdFh0Vnz2Q+4oa1RKdKJ4cQIRwMXdJBYWGQF8NMXL2TH50WQ04MS7GYAYsOF0OYlDvtYliLOmD1kDbCROfKlrIzd0U5KANChjVFh4jFQmec78gj8oEPeZAOnXsM/8XwJ2KHsGLpkB9bHEOCWEycY4iWeTQcc319+vZxkUFs4YZlF8OgHij5Q8eOmMOJuiMAsOACAYZ08uHYp/6chwFWK5YQ3HHkiR9huEiZklh448ePd7HnHHXiouTss8+uGa7t1NEQvkMPPdTrNmbMGLdmGf6DWTj2ERW+K4g89eXChAsi0tpvv/3c4kOsuQAhPVgelrQz+TNUmO+COwI4bNgw/47QVrvuuquXLcQzP56OS5OABKw0263FS01HQMeOdUAHS0dK50sHiajQoTLcxFU9HTeChXUWV8Gcp1OhY8UKohNFyKKDo9OlE6XDpIPkyp/ho2eeecY7MYbIEBT8iEv+dICxLJ1OGKuNeORFeQuJTYDEj/mh+FAmHB03QsEiAFYtPv/8815u8goXnWYcp/MJPwSW83T44Q8XeCD8WFUs1oAlnT9zZYgb5SEc5aGOCOKG9Rtc7DkmLT6EDccxggBDHGUdlAxnYiGFZRhhY0u6fKL+pEeeUV7OUxfOMU82OBkSxPLhw4UH7YsQE4bFK7QtYsOHNsXqpB5p53VJ8mEb5aK9+A7FUB8XFKSPRccFDOUgbb5bOOoajNPHrC7lAosLGi46XnvtNb844Tsh13oISMBaT1s2a03oOJiDYtiNzjLfRedGODoNLI7lK5a74BCWjggrKSwRwtEx0tHj6NAYLqKzYuiK+QzCf+tb33JRuemmm3yBBPFwdLhc4bMlbzpM9jt1Tib/NzeSanV6nkDyh44wOkPSwCF+++67r88RYYVhIWAJUHdchPeDTX+i4w5/RBhrEU7EYw6JeJQdoedCALGmfszXMEeIkBEGkWKLSLONfY7Dxfn0cZolIo94IcbUh7hpMYFvCHakxTbd2XOMIy7pUNZgT16ky3naLcQuysMWkcp36bxghiMdykNb4mhH0ufCJMSUMkR4D7TpD+UhTerL8DKWMkOYDH8yxMl3aVAi5BEuHVf7pUmg5ltTmmVXqVuYAB0VHR+dE51CupOKzo0i0qlgBfTt09cXRDAPFMunR48e7ZYC6SAOLGigg8fS4phhIDouFo1w1Y1gYunRyUUnHB0S4adPn+5X3Oxz9b7r4BpLgLLiolyx9ZObztMxMvSEFYB4RvpYBCx44D41hrOwFBHe/DRIK5EYj0f5sEgY3kNsmRdiTgsOxIMVIonlQt3ooBkSY+iUeTE6cDpprCgsWm7URpCJG+IaZWeb36Fff/31XlbEk6FYLFasPdghosx50cEjyiyUQHxIG9aICKKxPGEQjvTJn+E+hjOZB5yVzDlyjxiLQBg6ZrgTV6h8nC/Ei/O48EO44cZiHMrO9yA9zBzhamLV/sv3DK60Eby5pYNynXTSSZ4W5aRdCbeldGqnqqMsE5AFluXWyXjZuNquKN/8apgOgs6ObTg6bqwnbijm6hiHGDHnRcdHx0Nnw4o/FiTQ4bIIgDkc5mW4R4jJexaGICh08mzDITx0piyCYMiIOSfuUyMMaSFGlCc6r7TYkgbHiAjzcyzgYCHJgQceaCwsoNPnPJ0iQ1Ks9KvLcfMzq/MQB4QARqTLjbrUIYQGQWUfobvrrru8XoSnnMz/MGyGGB5xxBG+cAUmLIiBFdZUvps1e5afon7UhXgsDkE8YXVOsgIP4UGoDkvmkFggw71otAGdOryIR92oK4sfHn/scTvjzDP8wgFxp03hiSCwmIb7rIjDEOGXvvQl3zJ0R53TjjC49PehkD/nmHc788wzDQHme0LZ+I5QH3iRVv53i3gIEumTP6szmcPk+wRTLC/qf8kllxgXTMyt1VUW0pIrHQJ6lFTptFUmShodBVYSV8iDBw92q4XChR9X0IgGFhAdXnTadNYsJGD4ieEpRC3mxLDK6Kzo2FeuWOn3FJE2c1nExwrhip80sEq4yudDfDojFhrQcdK5MYSE6BAfS4kOn04cy4OOmryw+BjSIz3mR7B8SJeyYx2SBqKItcUwH3MwdIKsgGQ1JEN7hRwiSl6w4EM40mWIkLxxlIcywIc8EFxYUnYEhrCwobOmQ6be+HGONGCLZRrpUQeGKAclw2OIPRZdsIApVhWdOPnhsO5IE+sPiwpxSrTSBg0c5GILa+rBnBOCRjzKS1vSHogs8VlsQR3xhzX50LbkT/kQW8pNm0dbESbtyIsPbUOb4UiD9GkL+GF5Uw54UHfqi0Uc36t0etSNssOVeF63JAAcwoqnvBKwNLXS3ZeAlW7bqeTNRIBOkRWU3NPGEn0smRDrhhRha+I0JP3WEFaMWkMrNl8d2swQ4mY/DCb2Y4QrvR/sN52LeL4lQl6c8M9FS65IwzEElAsfJ9kWyo/TSdz0lWH6OPbTW5Ly8El6yfV+rbj41XLpPNP7mwKl002XoVYayUE6XOwTPvYJz34cux9lS/5xLvzZRj4RP+JEOPd3hDVxa533lDb/kx8ml/amMqRjJMVOyvRRnXJhN5WTK3zmgriJl6HJX/7yl37zcjoN9qPcsR9liHBej+TAvw+bwodfbNNxCB9lwT/8Ip3IL32cPhfh+U7wHxdpsp/z31TPSCfC4R9hwo9jT4N23JRofhj8I66H3XRMuuEivTiONDgOv0gjwkRaHBcKH+Hy/TnOTzOdVjpe/n7k4/FhGL/7/ICFjrcQnnTr7BeStNyfL2Z9LpXHluIU8qt1LpWOZ5k+Tu9vKk+tuFHGAuHCq6m3ssCamrDSL2kCrIJjIQHDZszPxDBcSVdKhReBVkKgTQgYVw2MnxcaM28l7ahqNBGBuIKP5Pke5Z8LP21FoC0SoH9tqb61TQwhTkmeHMANs0zG0wHJiUBDCMR3ptjhp4akrbAiUMoE+E2wUIdbFlgIVXCIsQkr2CYEjGfAsRIKyMDWFXQTfqOUtAiIQKsnEFYXDx74+c9/7o8tQ8Ca27UJAWPpLs+G4yGlciIgAiIgAo1DgLlhnlNZ183rjZNL3am0CQHjagHLC8e+nAiIgAiIwLYRYCQLAePTUq5NCFgaroYP0zS0LwIiIAINJ5AVQ0DPQmx42ymGCIiACIhABghkygJLq3pdllI6DPzqCpcBtiqCCIiACIhAExLIlIAVI0bFhGlCXkpaBERABEQgIwRaXMCwqBAlHp7Ke3t4QChPt+Ylgty3Ff7w4kGvLInnYavE4cnSAwcOzD2wMyNMVQwREAEREIFmINDiAhZ15JE9vCqD9zghZry6AgELh5Dx1O5rr702dy8XT+D+whe+4E8kTwsdcTjGIXSx7yf0RwREoFkJ6PfXrLi3OrNSHN1qcQELaLx64ayzzrIxY8b4e5+iFUKYeP0FL7njfoPzzz/fH13y7W9/21+dEK/kiDhsI13202+W5VhOBESg+Qikf4vNl6tyagsEWlzAAjIvnkOIeK9TofsKELDXX3/d31rLe4MQsr59+/rbc3n3Uf77mRiK5P0/WHbcaMdL+uREQASalwC/5WnTnk2mBuZbOaMh8Xj85i2GcquDAE/Gr0pGq3onb0vf92P7Jv1ouzpCZvN0ZgQs8NQ13MB5hhZ5uy3ixXG3rt1coELAwlojLV5gx3wZL8vjhYXjxo2LLLQVARFoYgLxW+Tlkz/92c/s//51d5Ijv8G1TZyzkm8YAV4w+oSNP+54u+H66xIB61lr3UHD0mr+0JkSML70fHD5W85hZXFFF0+WX7e+xvKKt66mhyp4LteQIUN8KJFXYMQDWUlHTgREoHkI8LtrX1lte53xE+s/eGQyf70uybiI9101T/HaeC7ViTHQ3ua9/Yp1WDXVLbFSA5IJAYurNYSI5xYyZxWvHg9x4hxDjLxSnFedY4XxWnvmzngNe6QRDYCYRdzYhp+2IiACzUUgGaJKVg+vqSq3DWWVllx+1pqfbq5SKJ/NCbiRkLTJmuqaNto8RPbPtLiAhfAwV8UcFws1ZsyY4asRR44c6QQRK+a9xo4da1dccYXdf//9/iNgOLHQAo587N5Q+Sd1LAIi0GwEmGex6uRVRoyybDFXfLNioTVFWbKQ5qYy+IhXlV/8l2of2eICFt/llStX2uTJk42l8bz5dtKkSb5cHusK64tFGCNGjLCjjz7a3+1VVl5mp512mg0ePNiTSA8fRpraioAIZIFAcq+nF6MYYSomTHPVqSnKkoU0PyoDizhK2bW4gIXwME81YcIE27BxQ3KhVu3Df7H8PYYD2Z544on2iU98wpkzJ4Z1JicCIiACItD2CGSm90ecWEpfn0PUQtjqCyt/ERABERCB1ksgMwIG4kLjsGGhpZsgwhXyS4fTvgiIgAiIQOslkCkBK1aQig3XeptNNRMBERABEdD7wPQdEAEREAERKEkCErCSbDYVWgREQAREQAKm74AIiIAIiEBJEpCAlWSzqdAiIAIiIAISMH0HREAEREAESpKABKwkm02FFgEREAERkIDpOyACIiACIlCSBCRgJdlsKrQIiIAIiIAETN8BERABERCBkiQgASvJZlOhRUAEREAEJGD6DoiACIiACJQkAQlYSTabCi0CIiACIiAB03dABERABESgJAlIwEqy2VRoERABERABCZi+AyIgAiIgAiVJQAJWks2mQouACIiACEjA9B0QAREQAREoSQISsJJsNhVaBERABERAAqbvgAiIgAiIQEkSkICVZLOp0CIgAiIgAhIwfQdEQAREQARKkoAErCSbTYUWAREQARGQgOk7IAIiIAIiUJIEJGAl2WwqtAiIgAiIgARM3wEREAEREIGSJCABK8lmU6FFQAREQAQkYPoOiIAIiIAIlCSByqyVurq6ulaRysrKah1zUEyYzSLphAiIgAiIQKsikDkBKyRY+cSLCZMfR8ciIAIiIAKti0CmBGzjxo32zjvv2CuvvGLsDxkyxIYOHWrt2rXLUV+1apW9+OKLtmjRIisvL/cwu+66q7Vv3z4XRjsiIAIiIAKtn0AmBIwhQayqxYsX23XXXWfz5s2ziooKe+ihh+yCCy6wQYMGeUusX7/ennnmGbvlllusc+fOLnJTpkyxr3zlK7bbbrv50GJYZzHMyHHst/7mVA1FQAREoO0QyMQijhCd2bNnu0B961vfsu9+97u2bt06t7bWrFnjLfLhhx/a1KlTbZdddrHvfOc7NmHCBPefP3/+Zi1GmpFu2oLbLKBOiIAIiIAIlCSBFrfAwvpavWq1zZ071y2pHXbYwYcNR44caa+//rodeOCB1qlTJz+3++672+OPP26TJ082hhP33ntv692792bw33//fXvrrbcM0Zs+fbrtu+++m4XRCREQAREQgdIlkAkLDHzrN6y3lStXWvfu3Z0m1lO3bt1cpJgPw2FJMS9GuLvvvtvuueceD9+lSxf3T/9ZsmSJPfX0UzZx4kR77rnnctZYOoz2RUAEREAESpdAi1tggY4FGZWVlcY8Fw7LbMOGDX4uwjCUiPU1duxY+8QnPuFCduGFF9rMmTNtwIABEcy3LP5gcQfuiiuusKqqKt/XHxEQAREQgdZBoMUtsJin6tixow8Fzpkzx4ULq4s5sX79+llYWIjbtGnTbMcdd/Tz+BFv9erV3hqRFgcIIhZbfFpHc6kWIiACIiACQSATFhjWFkIzKFltiKXE6kOWxb/55pt25JFH+pJ5rC/CjBo1yhd6MLzIOea4Yg4s5tOicrHlvJwIiIAIiEDrIpAJAQukWFSf+9zn7I477vAViAwTDh482F544QXr2rWrjRkzxk444QRfRn/TTTe5oJ155pl+r1ikoa0IiIAIiEDbIJAJAYuhP+7tOvzww23cuHFOHysMqwsrjDDMkSFo559/vt8DRqAOHTrkbmKOdNpG06mWIiACItC2CWRCwNJNUGjOKv2UDUQKoZMTAREQARFo2wQyJ2A0R3rOqpBVVZ9/227Stlf79Peh7dU+mzWmTfjtVldXJZ9sllGlKn0CmRSwQqKVRl2ffzqs9ls/AX0fstfG0Sbt23fwFcHSsOy1UWsoUSYFrDWAVR2aj8DSpUuNJ68ktruu9psPez05VSfPM600HiiwfPkya985scbqiSFvEWgoAQlYQ4kpfIYI0CWW2f33T7TPfvZ06zNgV1u8fFVyK4a6yiw0UjKImCy8qrCqpfNsxNkn2EaNJWahWVpVGSRgrao521plal52umzp4qTih9iwE79gazfoiSvZ+BYk8lVeYevXrLSZ9/4uKRJyxl85EWg8AhKwxmOplFqIgMvYwB2s43Z9rWJjMoyYlKNG2moXKN155vun40S4CJP2q51izVF++EJhOFdXOoXixzniUY503NiPLWHSLn0+vZ8Ow36+H8e4YutdE7rwX1/EkbwSqbLD8uSTPKtU1ldhUDq7TQQkYNuET5GzQMA73rUbrCp5duZGnnkZPXEWCtdGy0ATlCeiVbVxg1XrOaRt9FvQ9NWWgDU9Y+XQHAQwG9x0SP6ECdEc+SqPggR8yUauHXI7BcPqpAhsLYEWf5jv1hZc8URABERABNo2AQlY225/1V4EREAESpaABKxkm04FFwEREIG2TUAC1rbbX7UXAREQgZIlIAEr2aZTwUVABESgbROQgLXt9lftRUAERKBkCUjASrbpVHAREAERaNsEJGBtu/1VexEQAREoWQISsJJtOhVcBERABNo2AQlY225/1V4EREAESpaABKxkm04FFwEREIG2TUAC1rbbX7UXAREQgZIlIAEr2aZTwUVABESgbROQgLXt9lftRUAERKBkCUjASrbpVHAREAERaNsEJGBtu/1VexEQAREoWQJNJmC8UlxOBERABERABJqKQJMJWFlZzVtYETKJWVM1n9IVAREQgbZLoEkE7M0337THHnvMli5daghZLTEzWWZt9+ummouACIhA4xGobLykzC0txOq9996zm2++2Z555hkbNGiQDRs2zHbddVfr2LGjZ4dFFqLWmPkrLREQAREQgbZDoFEFLERpzz33tFNPPdVeeeUVe/HFF33br18/Gzx4sIsZ+3IiIAIiIAIisC0EGlXAoiDbbbedHXzwwTZ8+HB7+OGH7corr7RHH33UTjnlFNtnn33skEMOsdGjR1tFRUVEyW3z58tCFHMBkp1iwqTDa18EREAERKD1EWgSAZs7d65Nnz7d2M6fP98+/vGP25e+9CXr06eP3X///XbttdfaXnvt5QKWP5xYSLDysRcTJj+OjkVABERABFoXgUYVsBCjd955xx544AEbOnSonXDCCW6Jde3a1ee9tt9+ex9SzFlfNYsVnerGjRtd8F577TWrqqry+TPmziorPyomeTDH9sYbb9iqVausS5cuLoa9evXKzcG1riZSbURABERABAoR+EgZCvk28FxYRojTJz/5STvyyCNzizWwxBYtWuRiw6KO8vKaBZDJGsWc8CxZssSuv/56YxUjaXTv3t0mTJhgu+yySy7M8uXL7c4777Rp06ZZj+16WI/uPWzAgAGGgMmJgAiIgAi0HQKNKmBYTVhRLN5AjA488EC3pBCjxx9/3J588kn72c9+Zu3atcsJEqhD+GbPnu3hfvWrX/mKxZ///Oe+CKR3795+vGHDBnvppZf8c8EFF9iIESNs9arV1q59O2+xSIeDmCfjXOx7IP0RAREQARFoFQQaVcDef/99n/tCrFasWGE77LBDTpxYjYillO9i2HH16tU+ZzZkyBCfK0PkRo0aZQwnHnDAAS5ga9as8aX5O+20ky1YsMA+/PBD69u3r3Gc79JiRlpyIiACIiACrYtAowrYunXrfA4LS4zhv3nz5jktxAQxQohiPistMARav369ix4rGHH4kwZCheWFI11EksUhixcv9jkwluSfffbZfp9ZiCFhsQCZi1u7dq2L4MiRIzktJwIiIAIi0EoINKqA9ejRwz796U+7WDHnxVxXDN8x54V/vnAFx7LEH3ELsSIew5G5xR6bAhK/sl2lffGLX7Ru3brZhRdeaC+//LLPk4U4EnThwoW+4vHdd9+1qVOnepkiL21FQAREQARKn0CjChiLL0JEeArHjBkzcnNdDBGOGTPGjjvuuFqiFILWKXlKx4477miTJk1yEUPw5syZ48ODrGDEkTbDkBwzdIggMj/GcCXWX+RNWMTz/PPP9/z/8pe/5IQUPzkREAEREIHSJ9CoAsaS9ligkV7+jkgxRMhKwRCsNDqsLeINSh47xZAfFlOHDh3s9ddft8MOO8yfqch8FxbX2LFj/TFVTz31lA8xYukdd/xx1qlTp3SSLpJhveEXlmCtQDoQAREQAREoWQKNKmBYUDzAF+tpv/32y0FBtGJerJCARcCYz7rtttvcojr22GNd1LgpGnEkTW6AxtK74447XBSPOeYYGzVylAsjIpVOP47JW04EREAERKB1EWhUAWP5PMLSuXPnWpQQNIYQx40b53NkYRlFoBCdzl0621FHHeXL7xEfHv6LJcbNz4Rp3769b5lnGz9+fPJc+yRMh44ejrQinfx041hbERABERCB1kOgUQWMR0UddNBBPheFSMWwHcLC4oydd955M5FJo+SmZkSKT9rFU+zjHEOC+UOG4aetCIiACIhA2yDQqALGECIP8cV98MEHtcQKMWNZfL6VVAhzCB9+hcKn/esKUyhdnRMBERABEWg9BBpVwBhCZDEGc05XXHFFTrAQIR4BdcQRR9hnP/vZWqsQC6EsJFrpcPX5p8NqXwREQAREoHUSaFQBY4iwvKzc56bOOusst55Y3h5PwmCIkfkwOREQAREQARHYVgKNKmDco4V1xA3IiBZL3JctW+aLOrhvi3u2ZD1ta5MpvgiIgAiIAAQaVcBidSEP5b3xxhvtiSee8Oca8lgnnij/uc99zl9oKRHTl08EREAERGBbCTSqgEVhXnjhBX+U009/+lPr2bOnP3T3nnvusfvuu894JiHWGQsxJGRBTFsREAEREIGGEmhUAWPoEFHiaRoMJ/I0Dpa7s6hj4MCB/hqUhhZQ4UVABERABESgEIFGFbBZs2b5UzJ4cSUP2P3jH/9ou+22mz9BnneBHXrooTmrS9ZXoebQOREQAREQgWIJNKqA8fDd22+/3Z+MwXzYzJkzfSiRlYh84kkaxRZO4URABERABESgLgKNKmB77723P20D8WI4EdGKYUWGEXnXl5bR19UUOi8CIiACItAQAo0qYDxtng+OB+6+9NJL/nBfjpkX22OPPax///713shMeDkREAEREAER2BKBRhWwWFnIfNddd91lb731lj8DkQUdDz/8sE2YMCH3fMQtFUp+IiACIiACIlAfgUYVsMiM+8B4d9fhhx9uvBGZtyePHj3ab2qOMNqKgAiIgAiIwLYQaJLnOjHfxWOjuHmZF1nyehXmv3gyB364/AfybkslFFcEREAERKDtEWgSC4whQx4hxbMRKysrfeiQ16mccsopufkvLaNve1821VgEREAEGpNAowpYiNKoUaN8BSI3MfP0+bffftstsOHDh+cErDErobREQAREQATaHoFGFbDA17VrV3vttdfssccesxUrVliPHj38Qb7Mi8mJgAiIgAiIQGMQaNQ5sJjX4lmIV111lc2bN8+fhbh48WL7+9//bpMmTfL7whqj4EpDBERABESgbRNoEguMe8C4Yfm8884z3tK8atUqu/XWW+3pp5+28ePH+zBiLLlv2/hVexEQAREQga0l0KgWWMyBseKQd3/xRI4OHTr40+e7d+9ugwYN0rMQt7alFE8EREAERKAWgUa1wB55+BF7YPIDfu/X/fffb9OmTfOnbyxcuND+/Oc/2yWXXJJbPi8LrFY76EAEREAERKCBBBpVwDp17uSWF5bWIYcc4vd88SxE7gcbM2aMv14lnoUY1loDy6vgIiACIiACIuAEGlXAxo4da3xwzHvxKCnextylSxcbPHiwbb/99u6nPyIgAiIgAiKwrQQaVcBiWJDVh/dPvN9eeP4Ff6TU6tWrbecBO9v448bb7rvvnpsH29bCK74IiIAIiEDbJdCoAhYYWUb/2NTH/Mkb/fr186dyMCc2ceJEGzp0qD+dI8JqKwIiIAIiIAJbQ6BJBIxhQ16dctxxx+XKhBU2ZcqUWs9C1DxYDo92REAEREAEGkigSQSsb9++NnXqVLvlllt8UcfKlSv9eMiQIbkXWkq8GthSCi4CIiACIlCLQKMKWIjSXnvtZQsWLHDRat++va1dt9aGDR1mRx55pJ6FWAu/DkRABERABLaWQKMKWBRi0aJFvuqQ94G99957vpBjp5128teqhMhFWG1FQAREQAREYGsINKqAxSpEHiU1d+5c22+//ax///5bUy7FEQEREAEREIEtEmhUAYucGDbkSfQ8hR4B4+Zl3gfG3NiIESNy82ARPr1FBNOuLostHa6uMOl0tC8CIiACItC6CDSJgPXq1cuHDSdPnuxi1bFjR1u+fLkdc8wxtueee25RwIoVo2LDta7mUm1EQAREQASCQKMKGKLC+794gC8vsuQeMJ6+UVVV5Z/OnTtv8R4wwjFnxhM82OfNzrzVmYcCpx3WF0v1Z86c6VZdoTDp8NoXAREQARFofQQaVcB4aO9DDz1kzz77rFtZ++yzjx1+2OHGAo4tuZg7++CDD+zGG2/0+Ighz1D86le/mptHi3Dr16/3m6LPPPNMf+/Yaaed5o+rSucRQ4ykE/tpf+2LgAiIgAiUNoFGfZ0K81533XWXjRw50ocKH3jgAZv84GR/iSUiglW1JTdnzhwXwO9///v2wx/+0K2xF1980dauXevRQoyw0Dh/wAEHWKdOnQoKFGH54Nq1a+db/REBERABEWg9BBrVAlu2bJkddthh9vnPf94Jde3a1RClGAIMQUnjC6uKJ3XMnj3bl98z9MhCkNGjR9urr77qqxkZlsQtXbrU7rnnHhs+fLhxg/TGqo05AYu0Itz8+fNt3bp1xqpIngwiJwIiIAIi0HoINKqAMbTHvBTPQmTlIZbS22+/7e8FQ1x69uxpPI2jkJARl/kzFoDgz4cXY77xxhu+ghHkiFGkPW7cOOOhweVl5QXn1fC788477Z133vEbqlk8IicCIiACItB6CDSqgCE+kyZNspdfftmH/bCasL4uuugiF5LzzjvPvva1rxUc0kOwED3eHxaOIUfOIX44LLwbbrjB3y3GzdJYbDie+jFw4EAP6yeSP8OGDbMLLrjA4/71r3/NpRH+2oqACIiACJQ2gUYVsOOPP94OS4YQER1cCA/ixD6Clj8fhR+OpfY77LCDz4EhYnywnvr06ZNboIGgdevWzXiy/YMPPmivv/66Pf300z7nxj1mzIeFq6yszFlmvI8syhL+2oqACIiACJQ2gUYVMJbJ82moQ1yY8xqUvMmZ+8WeeeYZF7sZM2YYVhsvx0TQGIK8+OKLfTEI537xi1/43BbDiWnxIv8QLASyvsUjDS2vwouACIiACLQ8gUYVsBCNLVUrLK5CYXhqxxlnnGHXXHONz3fxOhbe5Dxt2jRjQcj+++/v82LExRJjmBB//PLdlvLJD6tjERABERCB0iPQqAK2taIR8RjqGz9+vB100EFuQXGMZYXlxbAkQ5AhksytfeELX/BhwlihWHr4VWIREAEREIGtJdCoAra1hYh4CBmClT8cWJeFxbMW5URABERABNomgUwJGE1QnfzjfzhEDasrrLQ4zzassUJ+6XDaFwEREAERaH0EMidgyR1gxv+0q0ug6jqfjqt9ERABERCB1kmgUR8l1ToRqVYiIAIiIAJZJCABy2KrqEwiIAIiIAL1EpCA1YtIAURABERABLJIQAKWxVZRmURABERABOolIAGrF5ECiIAIiIAIZJGABCyLraIyiYAIiIAI1EtAAlYvIgUQAREQARHIIgEJWBZbRWUSAREQARGol4AErF5ECiACIiACIpBFAhKwLLaKyiQCIiACIlAvAQlYvYgUQAREQAREIIsEJGBZbBWVSQREQAREoF4CErB6ESmACIiACIhAFglIwLLYKiqTCIiACIhAvQQkYPUiUgAREAEREIEsEpCAZbFVVCYREAEREIF6CUjA6kWkACIgAiIgAlkkIAHLYquoTCIgAiIgAvUSkIDVi0gBREAEREAEskhAApbFVlGZtoJA2VbEURQREIFSJlBZyoVvrrJXV1c3V1bKp4EEysoQLrVPA7EpuAi0CgISsCKasaaTLCKggrQIgYqK5Gu8ERGTmLVIAyhTEWghAhKwOsBjdSFcixcvtqmPPWYfrl5VV8jkfHr4KqyB9Lk6om7V6fyOOvKJ83UlWp9/xCs2XISP7dbGi/gN3ZJfQr683J58/FGzgd1sY1VVQxNReBEQgRImIAGrp/Hmzp1rJ37600moscmnS/LZkHxCNJJduRYmUJHk/4gNOfzzNUZYC5dG2YuACDQfAQlYHazDAqvyq/r+dsQ3Jlh5p+2sukoCVgeyFjmNBTb76V1t7do1GkFskRZQpiLQcgQkYEWxX2brqsuswsqtqjpZuMnCgZoRrI9ih1FW33nCFROGlAmXDp/exx+XPhf7bHGRT6HzEaZQOI+c/MmPF2HjPMeRTrqs4Z/vR/hwET6O2ca52KbPxT7bSJ/dpD02RiGIJycCItBmCEjA6m1qesWNSedanfxP5ljqWpFYV+eZfz7/mPzzz6WP69qPchfyT59Lp58+n94vlFZd8dLn69qPtGObH45jXNq/5sxH59J+6f1a8RKP5MKiYDqRnrYiIAKtlkDmBCx/yXr+CsB8f1omP0xTtVZZcqVPfymXFQJqjKy0hMohAi1BIHMCVp8Y1effEhCVpwiIgAiIQPMTyJSAsWCCZeuzZ89ORuqqrW/fvtavXz+rqGClmRn+CxYs8M+GDRusffv2tssuu1ivXr2an5xyFAEREAERaFECmRCwWPG3dOlSu+mmm+zRRx/1YcHhw4fbl7/8ZRcxKK1fv96efPJJmzRpkm3YuMGWL1tuJ554oh133HG23XbJCsFN924RNoYasdhin/NyIiACIiACrYNAJgQsUM6ZM8cmTpxoP/nJT6xDhw72y1/+0l566SXbfvvt/bg8WTK911572f77729dunRxMfvnP/9pe+yxh40ePTqS8W16qLFdu3a1/HQgAiIgAiJQ+gRaXMDCalqzZo0PHQ4aNMh23nlnQ3QQpRkzZti+++7rAlZZWWm77bZbbkiRoUPib9yYrBLMc8uXL7f33nvP1q1bZ7NmzbKhQ4fmhdChCIiACIhAKRPIzNPoERpEB1HCesLa6tmzpy1btsyY7woX82H+iKepU12Y+vTpE965LU/QuPHGG+2yyy6zW2+9tdlWKuYKoB0REAEREIEmJdDiFljUDtHiU/Pki5qzWFcIWb5buXKlPfDAA/b888/b+eefbwMGDPAg6WHDYcOG2QUXXOAW2pVXXql5sHyIOhYBERCBEiewuTo0d4U23crDnBdzXfPnz3cRY1hw3rx51rt3b+vcubOXCoFiqHFqYnlNmTLFvvKVr9jYsTyj8KNFG36Q/GEIsnv37tajRw/r1q2bBCzAaCsCIiACrYRAiwuY3xycWFoI2KBk/mvJkiVuWU2fPt1efvllX6CBaK1YsSJ53t1amzZtmv3ud78zLCyGGBE5/NPWF22D9cYHl7bq/IT+iIAIiIAIlDyBzAwhQpKhwNNOO80uv/xyX3wxfvx4Gzx4sD311FO+TJ5l9U8//bQtXLjQt4gc1hVDhSzSQLBCyGJb8i2kCoiACIiACBQkkAkBC7Hp2rWrffKTn7SDDz7YxQhxYviQe7yYC+vYsaOdc845dvrpp+eEikUdWGK4SKdgTXVSBERABESgVRHIhIAFUQQIwYo5rzjPXFY45snkREAEREAERCBTAkZzxLxVNA2ilh4azPcnnCyvoKWtCIiACLQdApkTsEJilD6X3m87zaSaioAIiIAI5BNo8VWI+QXSsQiIgAiIgAgUQ0ACVgwlhREBERABEcgcAQlY5ppEBRIBERABESiGgASsGEoKIwIiIAIikDkCErDMNYkKJAIiIAIiUAwBCVgxlBRGBERABEQgcwQkYJlrEhVIBERABESgGAISsGIoKYwIiIAIiEDmCEjAMtckKpAIiIAIiEAxBCRgAyLqHgAAEAxJREFUxVBSGBEQAREQgcwRkIBlrklUIBEQAREQgWIISMCKoaQwIiACIiACmSMgActck6hAIiACIiACxRCQgBVDSWFEQAREQAQyR0AClrkmUYFEQAREQASKISABK4aSwoiACIiACGSOgAQsc02iAomACIiACBRDQAJWDCWFEQEREAERyBwBCVjmmkQFEgEREAERKIaABKwYSgojAiIgAiKQOQISsMw1iQokAiIgAiJQDAEJWDGUFEYEREAERCBzBCRgmWsSFUgEREAERKAYAhKwYigpjAiIgAiIQOYISMAy1yQqkAiIgAiIQDEEJGDFUFIYERABERCBzBGQgGWuSVQgERABERCBYghIwIqhpDAiIAIiIAKZI1CZtRJVV1fXKlJZWVmtYw6KCbNZJJ0QAREQARFoVQQyJ2CFBCufeDFh8uPoWAREQAREoHURyJSAYVktXbrU5s2b51bWDjvsYL1797by8o9GOquqqtz/gw8+sIqKCuvfv7/16NGjdbWKaiMCIiACIlAvgUwIGMKFVbVs2TK79dZb7b777nMB23fffe3cc8+1Pn36+DFhELe//vWv9vzzz7uwnXnmmXbsscda9+7da1U2hhmJw34cV1VXWVnyrz4X4WtGNCN8MrzJCGdZ7WHO+tLKnL/XIXOlaliBogm8LeKAJJL91lC/htHIXmjaIJolfj5eyjiZvSIXXaJW8/1K+sVcIyXNVVXTNtH3bYkHYfL71i2Fbyq/TAhYVG7unLn2r3/9yy666CLr0LGjXfqb39jLL79sPXv2tPbt29uaNWvskUcesQ0bNtiPfvQjW758uV199dU2dOhQGz16dE7kSA+44dq1a2d8cOVlH1lz4V9oi3WHq4nXPrH2KpNPuxrxS6VdKG7Wz1X5ly9hVISQZ7su/IgqrDJpm7KkvWijGle6PYx3KUnxy0v8O5b8GL1NqpI26VBZ0za0TxlXhCVeN34/uNbQRpVJn8bvh/6OPtbrlRrx8hMF/kT/SpzoWwsEa/JT8Ytv8ozqyiCU/MMPP7RZs2fZgAEDbMiQIQ5l7733tldeecXFafvtt7d169a55YVltvvuu9vatWtd1JYsWbJZ8itXrbQPlnxg69evd6utQ4cONmvWLNu4ceNmYQudoFw0zuwkjtlqW/7BfKv8cLVVVRUXv1CaWTiHrFckf6qS7aYLriwUa6vLUJb82BYvT9p53Rrr+UHN0PNWJ5aBiOVJ23CJtTHpI2u6yQwUaiuLQNusX7PSFi5ZYd2WvWcVS961qo3rk9Q+urjcyqRbNBq/H2qwodQbKPmGtatsZ0uXvW/LVq6wN9+cab169XIDIQRqS6CZ2lm9erXNnz8/6RfpUZrftbiARZURJ+a/dtxxR7eegMMc2MyZM12ECIf4LFq0yBAzxIXjbt26OcR8s3fW27N8OPLNN9/0dPF/9dVXc0OJkW99WxrmtNNOsrK1L5it42tbut9aLK4NCbO3Z71tvXfsbdtt16PFvnj1cS/ev8wGD8Ja7mrly6eWbvskXyu+88tXLLf5CxbaroMGWWW7ygZ/X4vn1hwhk99L8rsbdsx+Vlb+rpUtXVC67bMJFx17zRy92eD+/Uq8fahUme3SbaOVdR1gv/vd77zv3VTVejewoF+dO3euYSC0hMuMgFH5AJIGwbm0qytMfjiGFSdMmOBXE3F1kF4Mkk6zbexXW0V5ha1YudIuueRnduqpn7GxiSW7NrlwyGdXajz4hpTuZUUNbTqCDslF2QvTp9u111xj/3XhhckFxnbJRdqGJEDt34Dap+UIVFSU2/XX35B03FV21llnt4ILwG1nSf/KNA+uufuSzAhY+/YdHMLChQv9SwGUBQsWuBXWqVMnh8M4LasS33//fR8+ZHiQhR+dO3d2fzqBAMgVQUtdFXhhMvqna2KxduzYKeG6o/VMhgvkskWAEQjmf3faaSfr0qVLtgqn0jgBVj3TP9FWci1LoMUFDMFBeDp27GCDBg1y0ZoxY4aLz4svvmhf+MIX3IpiAQeThR/72Mds2rRpPk+GePFlYkixkCNduRoCsMACZQEMHGMujx9iiL5YtRyBaB/ag/ahnXBqn5Zrk0I581uJD/7qY2ootVQf0uIClv6S7LzzznbyySfbpZde6vNeRx11lA0ePNieeOIJ69q1q+2333526KGH2htvvGE/+clP/IeOwO26666eTD7E/ON0Xm1tP1gwd3j44Ye7JQsDzodfW2OSxfpyMXb00UfnRg/a9rB3FlvIfAFZLD/Xb6dl2yhZ1ZodM4WirEzmaFiowT5zACzSWLFihVsP3OvFF4YhRMKxjxmPuMkVRwCurBxq6eWvxZW27YViWJwFTQyLq3PMVvvz26FNGA3CxdRGtkrZtkqTKQFrW+hVWxEQAREQgW0hkKkhRCrCVU7accUT5+KKNI4jXJyPY23rJsCtB6tWrfJbELDCuIrUMFXdvJrLh3bBMsYCYw6MUQV9r5uLfv350Odwryr3nrJfWVnpVnI88KD+FBSiKQhkTsAK/Wjzz+UfNwWY1pgmneRzzz1nt9xyi82ZM8fnFE899VRj7pEfpbg2f6sH93fffdeuvPJKe+GFF2zcuHF+CwgiFv7NXzLlmCbARd+9997rTwqiTZh3P/HEE/0hC+lw2m9eAsU9V6l5y6TcGpkAPzgcc4k33XSTL4z5+te/7othWNHJlaXEq5GhNzA5bvnYZ599/KKCG+5jFWIDk1HwJiLAatB+/frZOeecY9/61rf89zJx4kTjoeJyLUcgcxZYy6FovTkjTojYO++847cpnH322cZjunhKCbcsjBkzxgYOHKir/Rb4CsSFA6sPWXWLkM2ePVsXFC3QFlvKknvyuIWHYXeG3P/973/76mhu5YmbeLcUX35NQ0AWWNNwzVyqXNG/9957/qyzWD3FFSU/QCwzuZYlwJxKtAtDvXLZIsBcV8fkBnPEi98LQ/AIl8SrZdtJAtay/Js1dzpGfohx1c8+llkMMTZrYZTZZgQYpgoX9xnFsbbZILBu3Xp76KGHjIcscE+q3kXYsu0iAWtZ/s2WO1eOLAqIFYhkzL10DFkxLCLX8gRoIy4u+BT72p+WL3XbKQEXGFOnPuoLOU466SQbO3Zs26l8RmsqActowzR2sbC2WG3I8AdP02bokHetMYyoZ7o1Nu2tS49h3hg+ZDm9XLYIPPPMM3bzzTfbgQceaAcddJBf/Gn0omXbqOL/Ja5li6Dcm5oAPzKu6rG2uI/lrrvuSq4kp7qQcSW52267aQFHUzdCHelH27DA5qqrrrIbb7zRO0nai4UDffv29WHfOqLrdDMRYP748ssvtz/84Q/eJrynkHkwHi7O04LkWoaAnsTRMtxbJFc6y8WLF/t70bDEeOL5sGHDvKOMjrRFCtaGMw3uLMdmXoVhXc4xrMuFBatDsZ7lWpYAvxdEi5fn0jY87ov3FfL74ZF3ci1DQALWMtyVqwiIgAiIwDYS0H1g2wiwFKNzhR+OoSq5bBBIt0uUSO0TJFp+q/Zp+TbIL4EssHwiOhYBERABESgJAlqFWBLNpEKKgAiIgAjkE5CA5RPRsQiIgAiIQEkQ0BxYSTSTCpklAtwM/uyzz9rMmTP99Sd9+vSx/fff3++pY54k5q1izoRjlsmzco2nmEeY9Jb61YrH1GQyVVnrnAf66LynnxeOIOl0P4rvPkmCH8XnjJwIlDIBWWCl3Hoqe7MT4GZjHiXEU/15MzhPZ+A1G9xbx71CIRgUjP04fv755+2BBx7w8sa59Db2c/ESpdnsHOmlzuOfPvbEN+Ub6Xx0blN5UvHDT1sRKFUCssBKteVU7hYhMH/+fLv99tvtsMMOs1NOOcVfavj000/7Da68iHLkyJF+o+suu+zi99stX77cn5f3yCOP+JP/eZzXqFGj/P4unoTy9ttvu1ANGTLE3wqwcOFC44kPvNySJ6Rg2fGIKd4Txk3oS5cuddHca6+9/J4+LDueZE84nstHfMRy0aJFfoMtT1AfMGCAvfXWWzZ9+nR/Ekv37t09L8ooJwKlTEACVsqtp7I3KwGsL4YNERReR9O5c2fPn31uOn744Yf9ZtcTTjjBEAeE5KWXXrJjjjnG3++FAHHDMjfF8mLR66+/3ocUeUIKT3RAfLDkuGEW8ZoyZYoPUQ4fPtxuuOEGt/gQJF7lgR83OVOm6667zi6++GIbMWKEW3lz5871+HFj9BFHHGEPPvig5zl69GgXQoZBcTHc6Af6IwIlRkACVmINpuK2HAE6e54hyZMX0g9A5lUoPPYJQeOVKPHkDCwy9nfffXd/7BBheJEoQ49/+9vfbOjQoX6MgFWUV9hzzz/nw5E/+tGP/AkPvDDx7rvv9rk1XtvBsyzPPfdce/311+2MM86w2267zdPGKkPUyOu+++7zuL169fJj4mOBIY5YaocccohbiPEU9fQwZcuRVc4isHUENAe2ddwUqw0SoLNHpHjcUzx0F1FjHoy3WrNFRDiHCz+28eE8jyPCAkLAEC8sumRqyhd5EB/rjXxY8MFQIFYWgjl48GAXSsQMa4tHGfGOKgSVMvEiTPLBymPIEBE79thjPT0eQEs6WHJYeYSVE4FSJyALrNRbUOVvNgIVlRUuKogFQ4l77LGHiw/WD3NZO/XbyYUNocBhraWH6rDUcFg/WGdYRbiwgvDnKfQIEOKDpcacGQLHJ+ITHuELxzEfhA0L8OSTT7aDDz7YvdesWeN+/fv3d2uNN3D/9re/9fQQRFY6Ip5yIlCKBCRgpdhqKnOLEGDFH/NORx11lK9EZKEF1s+jjz7qc2CHJQs7Jk+e7B+WzPPEf0QHaw3ReuONN+z+++/3IT0Wctx5550+V4UYMTzIEB/vmLrllls8H+bJDj/8cLfGGCaMV6xg6S1YsMAtPkAgUogmaTBfxnwXi0ewxsKSe/fdd10QEUIe4sx5dxKvGg76W5IE9DqVkmw2FbqlCMQQH+LAAg2sL4TjU5/6lK8uxLJiEQVP/cfqQZB4Yjkr/xC1V1991c+zmIKwLPRgZSNCyApGxIU0SYOViccff7yLH3GxmHh/G4KIALF4hHk1hJTz5DNo0CC37ChbLOagHG+++aY99thjLnxjxozxuTDKJCcCpUxAz0Is5dZT2VuMACKCJYRj+A5hY8s5/Nj3G5HLa/wQPM7jYkgw0iBsnItwbDkfw4aRJuEiTDrPSJct1hhhcIQhTuTFOY45LycCpU5AAlbqLajyi4AIiEAbJaBViG204VXtbSeAlROfQqmFFRR+HtZXTcSZmpWK+eHw3VK6H8Wue69Q/ELn6k5BPiKQfQKywLLfRiqhCIiACIhAAQKywApA0SkREAEREIHsE5CAZb+NVEIREAEREIECBCRgBaDolAiIgAiIQPYJSMCy30YqoQiIgAiIQAECErACUHRKBERABEQg+wQkYNlvI5VQBERABESgAAEJWAEoOiUCIiACIpB9AhKw7LeRSigCIiACIlCAwP8HpQNP4vuVyHkAAAAASUVORK5CYII=" } }, "cell_type": "markdown", "metadata": {}, "source": [ "### TODO, Part (a): Calculating the Cumulative Distribution Function\n", "\n", "Complete the following to calculate $F_X$ for the CDF of the probability distribution $f_X$, and then \n", "demonstrate your code by calculating and displaying the CDF for the random variable of $X$ \n", "from the Example Problem at the top of this notebook. \n", "\n", "You must print out the CDF and then display it using `draw_distribution`. You should get:\n", "\n", "![Screen%20Shot%202021-03-12%20at%2011.00.35%20AM.png](attachment:Screen%20Shot%202021-03-12%20at%2011.00.35%20AM.png)\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def CDF(fx):\n", " pass # your code here, return the list of probabilities for the CDF\n", "\n", "\n", "# your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part (b): Generating random variates by inverting the CDF\n", "\n", "The basic idea here is that the CDF is a function from outcomes to probabilities:\n", "\n", "$$F_X\\,:\\,R_X\\rightarrow [0..1)$$\n", "\n", "if we invert this function, we get\n", "a function from the interval $[0..1)$ into the outcomes:\n", "\n", "$$F^{-1}_X\\, : \\,[0..1)\\rightarrow R_X$$ \n", "\n", "We used the same \"pinwheel\" approach in a previous homework: just generate a random value $a$ in the range $[0..1)$ and look for the first bin which is greater than $a$; output the corresponding outcome.\n", "\n", "The next cell contains a test of the CDF function, and demonstration of this idea: uncomment the code and run the cell a few times to see where the random value ends up: the first bin that the red line intersects going left to right indicates the outcome that is output. Be sure you understand how this process works, and what number would be output, for each test. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "'''\n", "\n", "Rx10b = [0,1,2,3]\n", "fx10b = [1/8,3/8,3/8,1/8]\n", "Fx10b = CDF(fx10b)\n", "plt.figure(figsize=(8, 4))\n", "plt.bar(Rx10b,fx10b,width=1.0,edgecolor='black')\n", "plt.ylabel(\"Probability\")\n", "plt.xlabel(\"Outcomes\")\n", "if (Rx10b[-1] - Rx10b[0] < 30):\n", " ticks = range(Rx10b[0],Rx10b[-1]+1)\n", " plt.xticks(ticks, ticks) \n", "plt.title(\"Probability Distribution Function\")\n", "plt.show()\n", "\n", "\n", "plt.figure(figsize=(8, 4))\n", "plt.bar(Rx10b,Fx10b,width=1.0,edgecolor='black')\n", "r = random()\n", "plt.plot([-0.5,3.5],[r,r],color=\"red\")\n", "plt.ylabel(\"Probability\")\n", "plt.xlabel(\"Outcomes\") \n", "plt.title(\"Inverting the CDF: r = \" + str(round4(r)))\n", "plt.show()\n", "\n", "'''\n", "print()" ] }, { "cell_type": "markdown", "metadata": { "scrolled": false }, "source": [ "### TODO for 10 (b):\n", "\n", "Complete the following code template to generate random variates for\n", "a given random variable (represented by Rx and fx), using the code from Part (a). \n", "\n", "Demonstrate your code by generating $10^5$ variates from the random variable $X$ from Problem Two and display it using show_distribution(...). \n", "\n", "Hint: Compare with the theoretical distribution from Problem Two. They should be very close! " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "seed(0)\n", "\n", "def rvs(Rx,fx):\n", " pass # your code here, generate the CDF, a random variate in [0..1) and \n", " # then the random variate required by the CDF\n", " \n", "# your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Eleven: Creating Random Variates for Standard Distributions by Inverting the CDF\n", "\n", "Now we will apply the technique from the last problem to generate random variates for two common distributions, which we studied this week:\n", "\n", " > Binomial B(N,p) \n", " \n", " > Geometric G(p): Recall that this is the number of flips you make until the first head appears with a coin whose probability of heads is p. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part (a) Binomial Variates\n", " \n", "Display the result of generating $10^5$ random variates with `rvs` just written in Part A from the Binomial\n", "Distribution with $N=8$ and $p=0.8$ using show_distribution. \n", "\n", "Hint: You may look at the Distributions Notebook on the class web site to see the\n", "formula for the PDF for the Binomial, or look at the lecture notes....\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Solution\n", "\n", "seed(0)\n", "\n", "N11a = 8\n", "p11a = 0.8\n", "\n", "# Your code here\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part (b) Geometric Variates \n", "\n", "Although $R_x$ is infinite, we will only approximate this distribution by considering the first 20 outcomes. \n", "This will potentially create an error, since of course it is possible for the value produce to be larger than 20, however, the probability is so small we will not worry about it. Note that our code for the CDF did not\n", "calculate the last bin, but simply set it to 1.0, which will make sure our generation of random variates does not crash. \n", "\n", "Display the experimental distribution for the Geometric Distribution with $p=0.6$, using show_distribution from the beginning of the notebook, for $10^5$ trials and the given N and p. \n", "\n", "Hint: You may look at the Distributions Notebook on the class web site to see the\n", "formula for the PDF for the Geometric, or wait until lecture....\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Solution\n", "\n", "limit = 20\n", "p11b = 0.6 \n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Twelve: Generating the Geometric Distribution by Explicit Formula\n", "\n", "Now we will explore using an explicit function for the inverse of the CDF. This is not possible for all distributions, but when it is, it the simplest (and most efficient) method. \n", "\n", "The following formula is from Wikipedia: if U is a random variable uniformly distributed in the range [0..1), then\n", "\n", "$$ 1 + \\lfloor \\,\\ln{( U )} \\, \\,/ \\, \\ln{( 1 - p )} \\,\\rfloor$$ \n", " \n", "is an integer which is distributed according to the Geometric Distribution with probability p.\n", "\n", "Note: $\\ln$ is log to the base $e$ (just log(...) in Python). \n", "\n", "\n", "For this problem, simply complete the following function stub and demonstrate it as in the previous problems. \n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "scrolled": false }, "outputs": [], "source": [ "def geom_rvs(p):\n", " return 0 # your code here\n", "\n", "# test it\n", " \n", "p12 = 0.4 \n", " \n", "num_trials = 10**5\n", "\n", "seed(0)\n", "\n", "# Your code here" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 2 }